1 Introduction

A binary constant weight code is an error-correcting code, defined such that all codewords have the same Hamming weight [1]. The codes have a wide range of applications, including the design of very large-scale integration systems [2], the frequency assignment problem [3], the design of optical orthogonal codes [4], error-tolerant key generation from noisy data [5], inter-cell interference mitigation for flash memory [6], and the design of Grassmannian constellations [7]. Furthermore, although not pointed out in the literature, binary constant weight codes have the exact same structure as index modulation [8] used in wireless communications and data storage studies. Using binary constant weight codes for codebook design is expected to achieve near-optimal performance due to their error-robust nature [9]. This codebook design is known as the index selection problem, which has been addressed in conventional studies [10,11,12,13,14,15] with some heuristic methods, and the problem itself is equivalent to the construction problem of binary constant weight codes .

A binary constant weight code \(\mathcal {C}\) is composed of M different codewords, each of which is a binary row vector \(\textbf{c} = [c_0 ~ \cdots ~ c_{n-1}]\) with a length n and weight w. Let \(d(\textbf{c}_m, \textbf{c}_{m'})\) be the Hamming distance between binary vectors \(\textbf{c}_m\) and \(\textbf{c}_{m'}\). Its minimum value \(d = \min _{0 \le m< m' < M} d(\textbf{c}_m, \textbf{c}_{m'})\) is termed the minimum distance of \(\mathcal {C}\). The maximum number of possible codewords is denoted by A(n, d, w), which is an upper bound on M. A concern in the coding theory [16,17,18,19,20] is to clarify A(n, d, w) or to construct a code that achieves A(n, d, w). It is widely known that this problem can be equally attributed to the maximum clique problem, which has been shown to be NP-complete [21]. Specifically, consider binary codes of length n and weight w as nodes in a graph, and connect each node with an edge only if the Hamming distance between these two adjacent nodes is greater than or equal to d. Then, the size of the maximum clique in this graph achieves A(n, d, w). Based on this structure, several studies have developed heuristics to construct better binary constant weight codes [22,23,24]. However, depending on the parameters (n, d, w), the optimization of binary constant weight codes becomes a challenging task due to its combinatorial explosion, even with the use of heuristics. Furthermore, proving that a code achieves A(n, d, w) generally requires an exhaustive search over the search space.

In 1965, Gordon Moore predicted that the density of transistors on integrated circuits would double approximately every two years, a historical trend known as Moore’s Law [25]. However, Moore’s Law is anticipated to be approaching its end due to physical limitations, and the computational performance of classical computers is expected to plateau. In response, quantum computing is expected to overcome current limitations and replace classical computing, at least for a selected number of tasks [26]. To offer some examples, quantum annealing (QA) [27] and quantum approximate optimization algorithm (QAOA) [28] are well-known heuristic approaches based on quantum mechanics. Although these methods have already been demonstrated to work effectively on real devices, both approaches have proven unable to outperform classical computation under realistic assumptions [29]. Consequently, there is growing interest in quantum algorithms designed for fault-tolerant quantum computers (FTQC) [30].

Quantum computation over FTQC has proven to have advantages over classical computation in terms of query complexity required for solving specific problems, such as Shor’s algorithm [31] and Grover’s algorithm [32]. Grover’s search algorithm finds a solution in an unordered database of N elements with query complexity \(O(\sqrt{N})\). Boyer et al.  extended Grover’s algorithm to the case where the number of desired solutions is initially unknown [33], which is termed the Boyer-Brassard-Høyer-Tapp (BBHT) algorithm. Then, Dürr and Høyer proposed Grover adaptive search (GAS) [34,35,36], an extension of the BBHT algorithm to solve optimization problems. GAS can be interpreted as a kind of quantum exhaustive search algorithm, which guarantees the optimality of the solutions obtained.

In the line of studies on Grover’s algorithm, an efficient construction method for a quantum circuit corresponding to any types of objective functions has long been unknown. The breakthrough is Gilliam’s method proposed in [37] that efficiently constructs a quantum circuit of GAS for arbitrary binary polynomial objective functions. This method has already demonstrated quadratic speedups of some problems, such as maximum likelihood detection [38, 39] and wireless channel assignment problem [40, 41]. In GAS, its stop policy and the Grover iterations can be appropriately tuned for a fast convergence to the optimal solution [42], demonstrating improvements in the average numbers of classical iterations and total Grover rotations.

Against this background, in this paper, we formulate the problem of finding a binary constant weight code as a polynomial binary optimization problem. This formulation enables the use of GAS and provides a quadratic speedup for the search problem. Additionally, focusing on the inherent structure of the problem, our approach further accelerates GAS, which is a unique attempt in the literature. Solving the problem with polynomial complexity is considered impossible [43], even with quantum computers, because the problem is equivalent to the NP-complete maximum clique problem. However, we believe that providing a quadratic speedup and reducing the often-ignored constant overhead of the GAS algorithm is meaningful as a first attempt. The major contributions of this paper are summarized as follows:

  1. 1.

    The problem of finding a binary constant weight code is formulated as a polynomial binary optimization problem, and the number of qubits required for constructing a quantum circuit is maximally reduced. In this formulation, the Hamming distance between codewords, which is usually represented by exclusive OR (XOR), is attributed to the representation using an inner product. Furthermore, a unique trick using the inner product is incorporated into the objective function to guarantee that the minimum distance of an obtained code is at least a given value d.

  2. 2.

    Two mathematical properties are clarified: (1) an upper bound on the minimum of the objective function value and (2) a lower bound on the exact number of solutions. These clarifications help reduce the number of qubits and query complexity, which are supported by our algebraic analysis and numerical simulations.

The remainder of this paper is organized as follows: In Sect. 2, we review conventional quantum search algorithms such as BBHT and GAS. In Sect. 3, we present our problem formulation, and our algorithm based on GAS. In Sect. 4, the proposed GAS is compared with the conventional GAS, where both rely on the proposed formulation. Finally, in Sect. 5, we conclude this paper.

Italicized symbols represent scalar values, and bold symbols represent vectors and matrices. We use zero-based indexing. Table 1 summarizes the important mathematical symbols used in this paper.

Table 1 List of important mathematical symbols
Full size table

2 Conventional quantum search algorithms

In this section, we introduce the conventional quantum search and optimization algorithms upon which our proposed algorithm is based.

2.1 Grover’s search and BBHT algorithm [33, 45]

Grover’s search algorithm [32] finds one of t solutions in N elements. The query complexity is \(O\left( \sqrt{N / t}\right) \), which is the total number of Grover operators applied to reach an optimal solution. Specifically, the probability of observing the desired solution is amplified by applying the Grover operator \(\textbf{G}\) to the uniform superposition state \(\frac{1}{\sqrt{N}}\sum _i {|{i}\rangle }\). When the Grover operator is applied L times, the success probability is expressed as [45]

$$\begin{aligned} P_{\textrm{success}}(L) = \sin ^2\left( (2L + 1)\arcsin \left( \sqrt{\frac{t}{N}}\right) \right) . \end{aligned}$$
(1)

According to (1), the number of Grover rotations that maximizes \(P_{\textrm{success}}(L)\) can be expressed as [45]

$$\begin{aligned} L_{\textrm{opt}} = \left\lfloor \frac{\pi }{4}\sqrt{\frac{N}{t}} \right\rfloor . \end{aligned}$$
(2)

Here, the number of solutions t must be known to obtain \(L_{\textrm{opt}}\), but in practice, it cannot be obtained in advance. To cope with this challenge, Boyer et al.  proposed the BBHT algorithm. As summarized in Algorithm 1, it searches for an appropriate L until the desired solution \(\textbf{x}\) is found. This repetition achieves the query complexity of \(O\left( \sqrt{ N / t}\right) \) even if the number of solutions is unknown. To be more specific, L is a random variable drawn from a uniform distribution [0, k), where the parameter k is increased by \(k_{i+1} = \min \left\{ \lambda k_i, \sqrt{2^{q_1}}\right\} \) for each iteration. The parameter \(\lambda \) is a constant value related to the increase rate of k. In [36], Baritompa et al.  demonstrated through numerical search that the expected value of the parameter \(\lambda \) that minimizes the query complexity is 1.34, while the optimal value depends on the distribution of the objective function values over the search space. They also found that when \(\lambda \) is set to 1.34, the expected total number of Grover operators is at most

$$\begin{aligned} 1.32 \sqrt{\frac{N}{t}}, \end{aligned}$$
(3)

provided that the \(\sqrt{N}\) ceiling of the parameter k, corresponding to line 5 in Algorithm 1, is ignored in the global optimization context, i.e., \(k \rightarrow \infty \).

Algorithm 1
figure a

BBHT Algorithm [33].

Full size image

2.2 GAS algorithm [34,35,36,37]

GAS is an algorithm that extends the BBHT algorithm to obtain the optimal solution for the polynomial binary optimization problems in the form

$$\begin{aligned} \begin{aligned} \min _{\textbf{x}} \quad&E(\textbf{x}) \\ \text {s.t.} \quad&x_i \in \{0, 1\} ~~~~ (i = 0, \cdots , q_1-1). \end{aligned} \end{aligned}$$
(4)

Here, \(\textbf{x} = (x_0, \cdots , x_{q_1 - 1})\) denotes binary variables and \(E(\textbf{x})\) denotes an arbitrary polynomial objective function. As with the BBHT algorithm, the query complexity is \(O\left( \sqrt{N / t}\right) \), where N is the search space size, i.e., \(N = 2^{q_1}\). As summarized in Algorithm 2, first, the initial value of the threshold \(y_0 = E(\textbf{x}_0)\) is determined by a random solution \(\textbf{x}_0\). Next, the BBHT algorithm is used to search for a better solution \(\textbf{x}\) such that the value of the objective function is less than the current threshold \(y_{i}\), and the threshold is updated by \(y_{i+1} = E(\textbf{x})\). Thereby, the optimal solution is obtained through multiple trials.

Algorithm 2
figure b

Conventional GAS [36, 37].

Full size image
Fig. 1
figure 1

Quantum circuit for GAS corresponding to an objective function \(E(\textbf{x}) - y = a_0x_0x_1 + a_1x_0 + a_2\), which is composed of \(\textbf{A}_{y}\), \(\textbf{D}\), \(\textbf{O}\), and \(\textbf{G} = \textbf{A}_{y}\textbf{D}\textbf{A}_{y}^\textrm{H}\textbf{O}\)

Full size image

Original GAS algorithm has assumed the existence of a black-box quantum oracle \(\textbf{O}\) that can identify the desired states. The breakthrough was Gilliam’s method [37]. In [37], Gilliam et al.  proposed an efficient construction method for a quantum circuit that corresponds to an objective function of binary optimization. Aided by two’s complement representation of the objective function value, the quantum oracle \(\textbf{O}\) can be easily constructed with only a single quantum gate. As an explicit example, Fig. 1 shows a quantum circuit corresponding to the objective function \(E(\textbf{x}) - y = a_0x_0x_1 + a_1x_0 + a_2\), which has quantum gates that correspond to the polynomial terms \(a_0x_0x_1\), \(a_1x_0\), and \(a_2\). Note that \(a_2\) is a constant that includes the constant term of \(E(\textbf{x})\) and \(-y\).

In the Gilliam’s construction method, \(q = q_{1} + q_{2}\) qubits are required, where \(q_{1}\) denotes the number of binary variables \((x_0, \cdots , x_{q_1})\), and \(q_{2}\) denotes the number of qubits for encoding the objective function \(E(\textbf{x})\). Here, \(q_{2}\) satisfies

$$\begin{aligned} -2^{q_{2} - 1} \le \min _{\textbf{x}} E(\textbf{x}) \le \max _{\textbf{x}} E(\textbf{x}) < 2^{q_{2} - 1} \end{aligned}$$
(5)

due to the two’s complement representation.

Next, we describe the procedure for constructing the quantum circuit for a binary optimization problem [37], where three types of unitary operators \(\textbf{A}_y\), \(\textbf{D}\), and \(\textbf{O}\) are constructed. The following steps 1–3 correspond to the construction of \(\textbf{A}_y\), while the step 4 corresponds to the Grover operator \(\textbf{G} = \textbf{A}_{y}\textbf{D}\textbf{A}_{y}^\textrm{H}\textbf{O}\).

  1. 1.

    Hadamard gates \(\textbf{H}^{\otimes q}\) act on the initial state \({|{0}\rangle }_{q}\) to create a uniform superposition state, yielding the transition of [37]

    $$\begin{aligned} \begin{aligned} {|{0}\rangle }_{q} \xrightarrow {\textbf{H}^{\otimes q}}&\frac{1}{\sqrt{2^{q}}} \sum _{i=0}^{2^{q}-1} {|{i}\rangle }_{q}\\ =&\frac{1}{\sqrt{2^{q_1 + q_2}}} \sum _{\textbf{x} \in \{0, 1\}^{q_1}} \sum _{i=0}^{2^{q_2}-1} {|{\textbf{x}}\rangle }_{q_1} {|{i}\rangle }_{q_2}. \end{aligned} \end{aligned}$$
    (6)

    Using the tensor product \(\otimes \), the Hadamard gates \(\textbf{H}^{\otimes q}\) are defined as

    $$\begin{aligned} \textbf{H}^{\otimes q} = \underbrace{\textbf{H} \otimes \textbf{H} \otimes \cdots \otimes \textbf{H}}_{q} \end{aligned}$$
    (7)

    and

    $$\begin{aligned} \textbf{H} = \frac{1}{\sqrt{2}}\left[ \begin{array}{ll} 1 & \quad 1 \\ 1 & -1. \end{array}\right] \end{aligned}$$
    (8)
  2. 2.

    Each term of \(E(\textbf{x}) - y\) corresponds to a unitary gate \(\textbf{U}_{G}(\theta )\). It acts on the lower \(q_2\) qubits to rotate the phase of quantum states. Let a be a coefficient of the term. Then, the phase shift is defined as \(\theta = \frac{2 \pi a}{2^{q_2}}\). That is, the phase advance represents addition, and the phase delay represents subtraction. The unitary gate is defined by

    $$\begin{aligned} \textbf{U}_{G}(\theta ) = \underbrace{\textbf{R}(2^{q_2 - 1}\theta ) \otimes \cdots \otimes \textbf{R}(2^{0}\theta )}_{q_2} \end{aligned}$$
    (9)

    and the phase gate

    $$\begin{aligned} \textbf{R}(\theta ) = \left[ \begin{array}{ll} 1 & \quad 0 \\ 0 & \quad e^{j\theta }. \end{array}\right] \end{aligned}$$
    (10)

    For a state \(\textbf{x}_0 \in \{0, 1\}^{q_1}\), \(\theta ' = 2 \pi \left( E(\textbf{x}_0) - y\right) / 2^{q_2}\) and \(\textbf{U}_{G}(\theta ')\) yield the transition of [37]

    $$\begin{aligned} \frac{1}{\sqrt{2^{q_2}}} \sum _{i=0}^{2^{q_2}-1} {|{i}\rangle }_{q_2} \xrightarrow {\textbf{U}_{G}\left( \theta '\right) } \frac{1}{\sqrt{2^{q_2}}}\sum _{i=0}^{2^{q_2} - 1}e^{ji\theta '} {|{i}\rangle }_{q_2}. \end{aligned}$$
    (11)

    The interaction between binary variables is expressed by controlled \(\textbf{U}_{G}(\theta )\). A specific method for interacting multiple binary variables is exemplified in Fig. 1, where \(\textbf{U}_{G}(\theta )\) gates are controlled by the corresponding qubits.

  3. 3.

    Applying the inverse quantum Fourier transform (IQFT) [46] to the lower \(q_{2}\) qubits yields the transition of [37]

    $$\begin{aligned} \frac{1}{\sqrt{2^{q_2}}}\sum _{i=0}^{2^{q_2} - 1}e^{ji\theta '} {|{i}\rangle }_{q_2} \xrightarrow {\textrm{IQFT}} \frac{1}{\sqrt{2^{q_2}}}{|{E(\textbf{x}_0) - y}\rangle }_{q_2}. \end{aligned}$$
    (12)

    In summary, the steps 1–3 are referred to as the state preparation operator \(\textbf{A}_{y}\). It encodes the value \(E(\textbf{x}) - y\) into \(q_2\) qubits for all the states. That is, it satisfies [37]

    $$\begin{aligned} \textbf{A}_{y} {|{0}\rangle }_{q} = \frac{1}{\sqrt{2^{q_{1}}}}\sum _{\textbf{x} \in \{0, 1\}^{q_1}} {|{\textbf{x}}\rangle }_{q_{1}}{|{E(\textbf{x})-y}\rangle }_{q_{2}}. \end{aligned}$$
    (13)
  4. 4.

    Finally, to amplify the states of interest, the Grover operator \(\textbf{G} = \textbf{A}_{y}\textbf{D}\textbf{A}_{y}^\textrm{H}\textbf{O}\) acts L times, where \(\textbf{D}\) denotes the Grover diffusion operator [32]. The oracle \(\textbf{O}\) identifies the states of interest that satisfy \(E(\textbf{x}) - y< 0 \Leftrightarrow E(\textbf{x}) < y\). Since we use the two’s complement representation, such states can be identified by a single Pauli-Z gate

    $$\begin{aligned} \textbf{Z} = \left[ \begin{array}{ll} 1 & \quad 0 \\ 0 & -1, \end{array}\right] \end{aligned}$$
    (14)

    which is acted on the most significant qubit of \(q_2\) qubits.

Note that an open-source implementation of the GAS algorithm is available from IBM Qiskit [47].

3 Proposed quantum search algorithm

In this section, we propose an objective function that attributes the problem of finding a binary constant weight code to a binary optimization problem. The problem is defined as finding a code such that the code length n, weight w, number of codewords M, and minimum distance d are each equal to the given values (n, w, M, d), which is hereinafter referred to as condition (n, w, M, d). Additionally, focusing on the inherent properties of the problem, we derive an upper bound on the minimum value of the objective function and a lower bound on the number of solutions. Using these bounds, we modify GAS and reduce the number of qubits and query complexity.

3.1 Proposed formulation

To obtain the optimal solution for a binary constant weight code A(n, d, w), we formulate the search problem as a polynomial binary optimization problem, which is supported by GAS. A challenging task here is that the objective function has to incorporate the Hamming distance between codewords. In general, the Hamming distance \(d(\textbf{c}_m, \textbf{c}_{m'})\) between codewords \(\textbf{c}_m\) and \(\textbf{c}_{m'}\) is expressed using the XOR operation \(\oplus \), which is not representable by a polynomial function. As long as the Hamming weight is constant, the XOR operation can be expressed by an inner product \(\langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle = \sum _{i=1}^nc_{m, i} c_{m', i}\), where \(c_{m, i}\) is the i-th element of \(\textbf{c}_m\).

Theorem 1

If a code has a constant weight w, then, for any pair of codewords \((\textbf{c}_{m}, \textbf{c}_{m'})\), we have

$$\begin{aligned} d(\textbf{c}_{m}, \textbf{c}_{m'}) = 2\left( w - \langle {\textbf{c}_m, \textbf{c}_m'}\rangle \right) . \end{aligned}$$
(15)

Proof

For any pair of i-th elements of codewords \((c_{m, i}, c_{m', i}) = (0, 0), (0, 1), (1, 0), (1, 1)\), the following equation holds:

$$\begin{aligned} c_{m, i}\oplus c_{m', i} = \left( c_{m, i} + c_{m', i}\right) - 2 c_{m, i}c_{m', i}. \end{aligned}$$
(16)

Since weights of codewords are constant, w, by summing (16) over index \(i=1,\cdots , n\), we have

$$\begin{aligned} \sum _{i=1}^n\left( c_{m, i}\oplus c_{m', i}\right) =&\, d(\textbf{c}_m, \textbf{c}_{m'}) \nonumber \\ =&\sum _{i=1}^n \left( c_{m, i} + c_{m', i}\right) - 2 \sum _{i=1}^nc_{m, i} c_{m', i} \nonumber \\ =&\, 2\left( w - \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle \right) \end{aligned}$$
(17)

\(\square \)

Theorem 1 allows us to transform the minimum distance of a code \(d(\mathcal {C})\) into

$$\begin{aligned} d(\mathcal {C})&= \min _{0 \le m< m'< M} d(\textbf{c}_m, \textbf{c}_{m'}) \nonumber \\&= \min _{0 \le m< m'< M} 2\left( w - \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle \right) \nonumber \\&= 2\left( w - \max _{0 \le m< m' < M} \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle \right) . \end{aligned}$$
(18)

Here, the following Theorem holds.

Theorem 2

The minimum distance of a code \(\mathcal {C}^*\) is at least d, where we have

$$\begin{aligned} \mathcal {C}^*&= \mathop {\mathrm {arg~min}}\limits _{C} s(\mathcal {C}), \end{aligned}$$
(19)
$$\begin{aligned} s(\mathcal {C})&= \sum _{0 \le m< m' < M} \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle ^l,~~\textrm{and} \end{aligned}$$
(20)
$$\begin{aligned} l&= \left\lfloor \frac{\log \left( {\begin{array}{c}M\\ 2\end{array}}\right) }{\log \left( 1 + \frac{2}{2w - d}\right) } + 1\right\rfloor , \end{aligned}$$
(21)

which is valid for the \(d \ne 2w\) case.Footnote 1

Proof

From (15), for any pair \((m, m')\) of a code such that the minimum distance is at least d, the following inequality holds:

$$\begin{aligned} \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle \le w - \frac{d}{2}. \end{aligned}$$
(22)

From this relationship, the upper bound on \(s(\mathcal {C})\) is given by

$$\begin{aligned} s(\mathcal {C}) \le \overline{s}(\mathcal {C}) = \left( {\begin{array}{c}M\\ 2\end{array}}\right) \left( w - \frac{d}{2}\right) ^l, \end{aligned}$$
(23)

where l is an arbitrary positive coefficient. Conversely, for a code such that the minimum distance is less than d, there is at least one pair \((m, m')\) that satisfies

$$\begin{aligned} \langle {\textbf{c}_m, \textbf{c}_{m'}}\rangle \ge w - \frac{d-2}{2} = w - \frac{d}{2} + 1. \end{aligned}$$
(24)

From this relationship, the lower bound on \(s(\mathcal {C})\) is given by

$$\begin{aligned} s(\mathcal {C}) \ge \underline{s}(\mathcal {C}) = \left( w - \frac{d}{2} + 1 \right) ^l. \end{aligned}$$
(25)

By comparing the obtained upper and lower bounds, if l is sufficiently large, the following inequality holds:

$$\begin{aligned} \left( {\begin{array}{c}M\\ 2\end{array}}\right) \left( w - \frac{d}{2}\right) ^l&< \left( w - \frac{d}{2} + 1\right) ^l. \end{aligned}$$
(26)

Thus, if l is sufficiently large, (19) is minimized, and the minimum distance of the obtained code becomes at least d. Taking a logarithm of both sides, the condition for l is given by

$$\begin{aligned} l > \frac{\log \left( {\begin{array}{c}M\\ 2\end{array}}\right) }{\log \left( 1 + \frac{2}{2w - d}\right) }. \end{aligned}$$
(27)

\(\square \)

Based on the above relationships, we formulate the problem of finding a binary constant weight code as a polynomial binary optimization problem, while maximally reducing the number of binary variables. We propose an objective function of

$$\begin{aligned} E(\textbf{x}) = f(\textbf{x}) + \rho \cdot g(\textbf{x}), \end{aligned}$$
(28)

where

$$\begin{aligned} f(\textbf{x}) = \sum _{0 \le r< r' < \left( {\begin{array}{c}n\\ w\end{array}}\right) } \langle {\textbf{p}_r, \textbf{p}_{r'}}\rangle ^l \cdot x_r x_{r'} \end{aligned}$$
(29)

and

$$\begin{aligned} g(\textbf{x}) = \left( \sum _{r=0}^{\left( {\begin{array}{c}n\\ w\end{array}}\right) -1}x_r - M\right) ^2. \end{aligned}$$
(30)

Here, for index \(r = 0, ~ \cdots , ~ \left( {\begin{array}{c}n\\ w\end{array}}\right) - 1\), \(\textbf{p}_r\) denotes the \((r+1)\)-th row vector of the combinatorial matrix \(\textbf{P}(n, w)\), which represents \(\left( {\begin{array}{c}n\\ w\end{array}}\right) \) combinations [44]

$$\begin{aligned} \textbf{P}(n, w) = \left[ \begin{array}{ll} \textbf{1} & \textbf{P}(n-1, w-1)\\ \textbf{0} & \quad \textbf{P}(n-1, w)\end{array}\right] \in \{0, 1\}^{\left( {\begin{array}{c}n\\ w\end{array}}\right) \times n}. \end{aligned}$$
(31)

In (31), \(\textbf{1}\) denotes a one column vector of length \(\left( {\begin{array}{c}n-1\\ w-1\end{array}}\right) \) and \(\textbf{0}\) denotes a zero column vector of length \(\left( {\begin{array}{c}n-1\\ w\end{array}}\right) \). In (28), the binary variable \(x_r\) indicates whether the \((r+1)\)-th row of the combinatorial matrix is used as a codeword or not. That is, the relationship between \(\mathcal {C}\) and \(x_r\) is expressed by \(\mathcal {C}= \left\{ \textbf{p}_r \big | x_{r} = 1\right\} \). The penalty function \(g(\textbf{x})\) ensures that only M from the \(\left( {\begin{array}{c}n\\ w\end{array}}\right) \) combinations is used and has priority over the cost function \(f(\textbf{x})\). Thus, the penalty coefficient \(\rho \) can be defined as

$$\begin{aligned} \rho = \overline{f} + 1 = \left( {\begin{array}{c}M\\ 2\end{array}}\right) (w - 1)^l + 1, \end{aligned}$$
(32)

where \(\overline{f}\) denotes the upper bound of f.

3.2 Reduction of qubits

Since the number of qubits required to construct a quantum circuit determines feasibility, it should be maximally reduced. Here, we modify the objective function (28) so that a single codeword is always selected, leading to a reduction in the number of qubits.

Theorem 3

If a code that satisfies the condition (n, w, M, d) exists, then for every possible codeword \(\textbf{c}\) in the combinatorial matrix, at least one code exists that satisfies the condition (n, w, M, d) such that the codeword \(\textbf{c}\) is included.

Proof

When we consider a code as a matrix where each row vector is a codeword, the Hamming distance between codewords does not change even if we swap any two columns of the matrix. Therefore, by swapping the columns of any code that satisfies the condition (n, w, M, d), we can construct a code that contains any specific codeword and satisfies the condition (n, w, M, d). \(\square \)

We define the matrix \(\textbf{P}'(n,w)\) by removing all row vectors from \(\textbf{P}(n,w)\) whose Hamming distance from the first row vector \(\textbf{p}_0\) is less than d. Then, the objective function with reduced binary variables can be defined as

$$\begin{aligned} E'(\textbf{x}) = f'(\textbf{x}) + \rho ' \cdot g'(\textbf{x}), \end{aligned}$$
(33)

where

$$\begin{aligned} f'(\textbf{x})&= \sum _{0 \le r< r' < q_{1}} \langle {\mathbf {p'}_r, \mathbf {p'}_{r'}}\rangle ^l \cdot x_r x_{r'}, \end{aligned}$$
(34)
$$\begin{aligned} g'(\textbf{x})&= \left( \sum _{r=0}^{q_{1}-1}x_r - (M - 1)\right) ^2, \end{aligned}$$
(35)

and \(\textbf{p}'_r\) is the \((r+1)\)-th row vector of \(\textbf{P}'(n, w)\). The number of binary variables \(q_{1}\) is equal to the number of rows in \(\textbf{P}'(n, w)\) and is expressed as

$$\begin{aligned} q_{1} = \sum _{i=0}^{w - \frac{d}{2}} \left( {\begin{array}{c}w\\ i\end{array}}\right) \left( {\begin{array}{c}n-w\\ w-i\end{array}}\right) . \end{aligned}$$
(36)

The penalty coefficient \(\rho '\) can be defined as

$$\begin{aligned} \rho ' = \overline{f}' + 1 = \left( {\begin{array}{c}q_{1}\\ 2\end{array}}\right) (w - 1)^l + 1 \end{aligned}$$
(37)

where \(\overline{f}'\) denotes the upper bound of \(f'\). Let \(E'_{\textrm{max}}\) be the maximum of the objective function value, i.e., \(E'_{\textrm{max}} = \max _{\textbf{x}}E'(\textbf{x})\). It is upper bounded by

$$\begin{aligned}&E'_{\textrm{max}} \le \overline{E}'_{\textrm{max}} = \overline{f}' + \rho '\overline{g}' \nonumber \\ =&{\left\{ \begin{array}{ll} \left( {\begin{array}{c}q_{1}\\ 2\end{array}}\right) (w-1)^l + \rho ' \left( q_{1} - M + 1\right) ^2 & (2(M-1)< q_{1}) \\ \left( {\begin{array}{c}q_{1}\\ 2\end{array}}\right) (w-1)^l + \rho ' (M-1)^2&(2(M-1) \ge q_{1}) \end{array}\right. } \nonumber \\ =&{\left\{ \begin{array}{ll} O\left( q_1^2w^l(q_{1} - M)^2\right) & (2(M-1) < q_{1}) \\ O\left( q_1^2w^lM^2 \right) & (2(M-1) \ge q_{1}). \end{array}\right. } \end{aligned}$$
(38)

Then, the number of qubits required to encode the objective function value is expressed by

$$\begin{aligned}&q_{2}' = \lceil \log _2(\overline{E}'_{\textrm{max}})\rceil + 1 \nonumber \\&\quad = {\left\{ \begin{array}{ll} O\left( \log (q_1) + l\log (w)\right) & (2(M-1) < q_{1}) \\ O\left( \log (q_1) + l\log (w) + \log (M)\right) & (2(M-1) \ge q_{1}). \end{array}\right. } \end{aligned}$$
(39)

Assuming \(w \ne 2d\), the coefficient l is upper bounded by

$$\begin{aligned} l&= \left\lfloor \frac{\log \left( {\begin{array}{c}M\\ 2\end{array}}\right) }{\log \left( 1 + \frac{2}{2w - d}\right) } + 1\right\rfloor \le \left\lfloor \frac{\log \left( {\begin{array}{c}M\\ 2\end{array}}\right) }{2} + 1\right\rfloor \nonumber \\&\quad = O(\log (M)) \end{aligned}$$
(40)

Using (40), \(q_2'\) of (39) can be simplified to

$$\begin{aligned} q_{2}' = O\left( \log (q_1) + \log (M)\log (w)\right) . \end{aligned}$$
(41)

The total number of qubits is

$$\begin{aligned} q_{1} + q_{2}' = O\left( q_{1} + \log (M)\log (w)\right) . \end{aligned}$$
(42)

The query complexity is expressed by

$$\begin{aligned} O\left( \sqrt{\frac{2^{q_{1}}}{t}}\right) , \end{aligned}$$
(43)

where t is the exact number of solutions and is analyzed later.

3.3 Further speedup for Grover adaptive search

3.3.1 Novel initial threshold

Let \(E'_{\textrm{min}}\) be the minimum of the objective function value, i.e., \(E'_{\textrm{min}} = \min _{\textbf{x}}E'(\textbf{x})\). Similar to (23), its upper bound can be derived as

$$\begin{aligned} E'_{\textrm{min}} < \overline{E}'_{\textrm{min}} + 1 = \left( {\begin{array}{c}M-1\\ 2\end{array}}\right) \left( w - \frac{d}{2}\right) ^l + 1. \end{aligned}$$
(44)

The objective function value \(E'(\textbf{x})\) may become greater than \(\overline{E}'_{\textrm{min}}\). By contrast, if it is smaller than \(\overline{E}'_{\textrm{min}}\), the obtained code has at least the minimum distance d. Therefore, the desired code can be obtained without updating the threshold of GAS. That is, we set the initial value of the threshold to

$$\begin{aligned} y_0 = \overline{E}'_{\textrm{min}} + 1. \end{aligned}$$
(45)

Additionally, if the objective function value is greater than \(\overline{E}'_{\textrm{min}}\), the obtained code is not good regarding the minimum distance. This fact helps to reduce the number of qubits \(q_2\). Specifically, the penalty coefficient \(\rho '\) of \(E'(\textbf{x})\) given in (33) is replaced by

$$\begin{aligned} \rho '' = \overline{E}'_{\textrm{min}} + 1. \end{aligned}$$
(46)

and we have an updated objective function

$$\begin{aligned} E''(\textbf{x}) = f'(\textbf{x}) + \rho '' \cdot g'(\textbf{x}) \end{aligned}$$
(47)

where \(f'\) and \(g'\) are defined by (34) and (35), respectively. Then, the objective function is upper bounded by

$$\begin{aligned}&E''_{\textrm{max}} \le \overline{E}''_{\textrm{max}} \nonumber \\&\quad ={\left\{ \begin{array}{ll} O\left( q_1^2w^l + M^2(w - d/2)^l(q_{1} - M)^2 \right) & (2(M-1) < q_{1}) \\ O\left( q_1^2w^l + M^4(w - d/2)^l \right) & (2(M-1) \ge q_{1}). \end{array}\right. } \end{aligned}$$
(48)

As with the \(E_{\textrm{max}}'\) case, the number of qubits required to encode the objective function value can be simplified to

$$\begin{aligned} q_{2}''&= \lceil \log _2(\overline{E}''_{\textrm{max}})\rceil + 1 \nonumber \\&= O\left( \log (q_1) + \log (M)\log (w)\right) . \end{aligned}$$
(49)

The number of qubits that can be reduced is bounded by

$$\begin{aligned}&\lceil \log _2( \overline{E}'(\textbf{x}) )\rceil - \lceil \log _2( \overline{E}''(\textbf{x}) )\rceil \nonumber \\&\quad> \log _2\left( \frac{\overline{E}'(\textbf{x})}{\overline{E}''(\textbf{x})}\right) - 2 = \log _2\left( \frac{\overline{f}' + \rho '\overline{g}' - 1}{\overline{f}' + \rho ''\overline{g}' - 1}\right) - 2 \nonumber \\&\quad > \log _2\left( \frac{\overline{f}' + \rho '\overline{g}'}{\overline{f}' + \rho ''\overline{g}'}\right) - 2 = \log _2\left( \frac{1 + \overline{g}'}{1 + \frac{\rho ''}{\rho '}\overline{g}'}\right) - 2 \nonumber \\&\quad = \log _2\left( \frac{1 + \overline{g}'}{1 + \frac{(M-1)(M-2)(2w - d)^l + 4}{q_1(q_1 - 1)(2w-2)^l + 4}\overline{g}'}\right) - 2, \end{aligned}$$
(50)

which is obviously positive since \(\rho ' > \rho ''\) when \(q_1 > M\) and \(d > 2\). Note that, for constructing a code that satisfies the condition (n, w, M, d), it is sufficient to obtain just one code whose objective function is smaller than \(\overline{E}'_{\textrm{min}}\).

3.3.2 Novel limit on the number of Grover rotations

The BBHT search and GAS algorithms try to find the appropriate number of Grover rotations, which increases the query complexity compared to the case where the number of solutions is previously known. Although, in general, the exact number of solutions t is unknown, the possible range of t may be obtained from the inherent structure of the problem. Our proposed algorithm exploits this information.

Theorem 4

If \(w \le d\) and \(A(n-1, d, w) < M-1\), then the exact number of solutions t is lower bounded by

$$\begin{aligned} t&\ge \underline{t} = \left\{ \begin{array}{ll} w! & (w-\frac{d}{2}=1) \\ \displaystyle {\min _{2 \le i \le w - d/2}} \ \left( {\begin{array}{c}w\\ i\end{array}}\right)&(w-\frac{d}{2} \ge 2). \end{array} \right. \end{aligned}$$
(51)

The proof is given in the Appendix A.

We improve the quantum search algorithm using the obtained lower bound \(\underline{t}\) of (51). Traditionally, the number of Grover rotations L is uniformly drawn from the semi-open interval [0, k), where k is a parameter that increases in each iteration of GAS. In the BBHT search algorithm [33], the maximum value of k is set to

$$\begin{aligned} \overline{k}_{\textrm{BBHT}} = \sqrt{N} = \sqrt{2^{q_1}}, \end{aligned}$$
(52)

which is defined on the basis of \(L_{\textrm{opt}} = \left\lfloor \pi / 4 \sqrt{N/t} \right\rfloor \) and the query complexity \(O(\sqrt{N})\). The same value is used in the original GAS implementation of IBM Qiskit [47].

Although \(\overline{k}_{\textrm{BBHT}} = \sqrt{N}\) is commonly used, its exact mathematical justification has not been provided, and the true optimal value is still unknown. In the following, we clarify the optimal maximum value of k, denoted by \(\overline{k}_{\textrm{opt}}\), that minimizes the query complexity. The probability of success \(P_k(t)\) is known as [33]

$$\begin{aligned} P_k(t) = \frac{1}{2} - \frac{\sin (4k\theta )}{4k\sin (2\theta )} \end{aligned}$$
(53)

where

$$\begin{aligned} \theta = \arcsin \left( \sqrt{\frac{t}{2^{q_{1}}}}\right) . \end{aligned}$$
(54)

Then, the optimal value \(k_{\textrm{opt}}\) that minimize query complexity is

$$\begin{aligned} k_{\textrm{opt}}&= \mathop {\mathrm {arg~min}}\limits _{k}\ \frac{k}{P_k(t)} \nonumber \\&= \mathop {\mathrm {arg~min}}\limits _{k}\ \frac{4k^2\sin (2 \theta )}{2k\sin (2\theta ) - \sin (4k\theta )} \end{aligned}$$
(55)

and is upper bounded by

$$\begin{aligned} k_{\textrm{opt}} \le \overline{k}_{\textrm{opt}} = \mathop {\mathrm {arg~min}}\limits _{k}\ \frac{k}{P_k(\underline{t})}. \end{aligned}$$
(56)

Theorem 5

If the ratio \(t/2^{q_1}\) is sufficiently small, then

$$\begin{aligned} \overline{k}_{\textrm{opt}} \in \left[ 1, \lceil \frac{1+\sqrt{2}}{2} \sqrt{\frac{2^{q_1}}{\underline{t}}}\rceil \right] . \end{aligned}$$
(57)

The proof is given in the Appendix B. In general, \(t/2^{q_1}\) is sufficiently small that it can be approximated as \((t/2^{q_1})^2 \sim 0\). The search range of \(\overline{k}_{\textrm{opt}}\) can be limited within the closed interval using Theorem 5. Then, we can search for \(\overline{k}_{\textrm{opt}}\) with negligible complexity; however, \(\overline{k}_{\textrm{opt}}\) itself is not expressed in a closed form. For example, within the closed interval, root-finding algorithm such as the bisection method [48] or Newton–Raphson method [49] is used for the first-order derivative of \(k / P_k(\underline{t})\), all the k that take the extreme values are listed, and \(\overline{k}_{\textrm{opt}}\) is one of the listed values or one of both ends of the closed interval. Note that, the lower bound obtained in Theorem 4 requires assumption \(w \le d\) and \(A(n-1, d, w) < M-1\). Thus, this extra technique cannot be applied to all possible parameters of the problem. Since there has not been much research on the number of solutions for constant weight codes, more general and stronger bounds may be found in the future. Moreover, Theorem 4 also holds for \(k_\textrm{opt}\) by replacing \(\underline{t}\) with t. This result suggests that the value of k need not be larger than \((1 + \sqrt{2})/2 \cdot \sqrt{2^{q_1}} \sim 1.207 \cdot \sqrt{2^{q_1}}\), even if the number of solutions is completely unknown.

The query complexity can be reduced by setting \(\overline{k}_{\textrm{opt}}\) as the upper bound on k. Note that, although this paper only focuses on the lower bound, an upper bound on t may be obtained for some problems, and a similar approach could further accelerate the quantum search algorithm more.

Overall, the proposed GAS is summarized in Algorithm 3. The conventional GAS is initiated with a random solution \(\textbf{x}_0\) and a random threshold \(y_0 = E''(\textbf{x}_0)\). In the proposed GAS, a strict threshold is set from the beginning, using the upper bound on the minimum of the objective function value, \(\overline{E}'_{\textrm{min}} + 1\) of (45). Additionally, instead of the conventional upper bound \(k \le \sqrt{2^{q_1}}\), we use the stricter bound, \(\overline{k}_{\textrm{opt}}\) of (56). These modifications further accelerate the efficient GAS algorithm. Note that, in [36], the probability of a success is maximized with \(\lambda = 1.44\) if the number of solutions is sufficiently smaller than the search space size (about 0.2\(\%\) in [36]). In general, if the initial threshold is as strict as \(\overline{E}'_{\textrm{min}}\), the number of solutions becomes much smaller than the search space size. Thus, we set \(\lambda = 1.44\) in Algorithm 3 following [36].

Algorithm 3
figure c

Proposed GAS.

Full size image

3.4 Number of quantum gates

We analyze the number of quantum gates that determines the size of a quantum circuit, which closely relates to the feasibility of the proposed algorithm. In the GAS circuit, the number of gates required for \(\textbf{A}_{y_i}\) has a dominant impact, which corresponds to the state preparation operator. Thus, we only focus on \(\textbf{A}_{y_i}\) later.

As described in the first step of Sect. 2.2, the Hadamard gate \(\textbf{H}\) acts once on every qubit. Thus, using (42), the number of \(\textbf{H}\) gates involved in (6) is

$$\begin{aligned} G_\text {H} = q_1 + q_2 = O\left( q_{1} + \log (M)\log (w)\right) . \end{aligned}$$
(58)

In the second step, the phase gate \(\textbf{R}\) is used to represent the coefficients in the objective function. Using (39), the number of \(\textbf{R}\) gates involved in (11) is

$$\begin{aligned} G_\text {R} = q_2 = O\left( \log (q_1) + \log (M)\log (w)\right) . \end{aligned}$$
(59)

The number of 1-controlled phase (1-CR) gates is the same as that of first-order terms in the objective function and is expressed as

$$\begin{aligned} G_\text {1-CR} = q_1 \cdot q_2 = O\left( q_1\left( \log (q_1) + \log (M)\log (w)\right) \right) . \end{aligned}$$
(60)

Similarly, the number of 2-controlled phase (2-CR) gates is the same as that of second-order terms and is expressed by

$$\begin{aligned} G_\text {2-CR} = \frac{q_1(q_1 - 1)}{2} \cdot q_2 = O\left( q_1^2\left( \log (q_1) + \log (M)\log (w)\right) \right) . \end{aligned}$$
(61)
Table 2 Number of quantum gates required by \(\textbf{A}_{y_i}\)
Full size table

Based on the above analysis, the number of gates required to implement the state preparation operator \(\textbf{A}_{y_i}\) is summarized in Table 2.

3.5 Specific example

Let us consider a specific example \((n, w, d, M) = (7, 3, 4, 7)\) of the proposed GAS using the objective function \(E''(\textbf{x})\) given in (47). In this case, the number of qubits is calculated as \(q = q_1 + q_2'' = 22 + 15 = 37\) from (36) and (49). The objective function is \(E''(\textbf{x}) = \textbf{x}^{\textrm{T}} \textbf{Q} \textbf{x} + 576\), where the matrix \(\textbf{Q} \in \mathbb {Z}^{q_1 \times q_1}\) is given in (65), and Fig. 2 shows the corresponding quantum circuit. Based on the combinatorial matrix \(\textbf{P}(7, 3)\), the reduced matrix \(\textbf{P}'(7, 3) \in \{0, 1\}^{q_1 \times n}\) is generated by removing row vectors whose Hamming distance to \(\textbf{p}_0 = [1\ 1\ 1\ 0\ 0\ 0\ 0]\) is less than \(d = 4\) as follows:

$$\begin{aligned}&\textbf{P}'(7, 3) \nonumber \\&\quad = \left[ \begin{array}{ll} 1~1~1~1~1~1~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0 \\ 0~0~0~0~0~0~1~1~1~1~1~1~0~0~0~0~0~0~0~0~0~0 \\ 0~0~0~0~0~0~0~0~0~0~0~0~1~1~1~1~1~1~0~0~0~0 \\ 1~1~1~0~0~0~1~1~1~0~0~0~1~1~1~0~0~0~1~1~1~0 \\ 1~0~0~1~1~0~1~0~0~1~1~0~1~0~0~1~1~0~1~1~0~1 \\ 0~1~0~1~0~1~0~1~0~1~0~1~0~1~0~1~0~1~1~0~1~1 \\ 0~0~1~0~1~1~0~0~1~0~1~1~0~0~1~0~1~1~0~1~1~1 \end{array}\right] ^{\textrm{T}}. \end{aligned}$$
(62)

The proposed GAS outputs one of the optimal solutions

$$\begin{aligned} \textbf{x}_{\textrm{opt}} = \left[ \begin{array}{ll} 1 ~ 0 ~ 0 ~ 0 ~ 0 ~ 1 ~ 0 ~ 1 ~ 0 ~ 0 ~ 1 ~ 0 ~ 0 ~ 0 ~ 1 ~ 1 ~ 0 ~0 ~ 0 ~ 0 ~ 0 ~ 0 \end{array}\right] ^{\textrm{T}}, \end{aligned}$$
(63)

which indicates that the optimal code consists of the first, sixth, eighth, 11th, 15th, and 16th rows of \(\textbf{P}'(7, 3)\) should be used, in addition to \(\textbf{p}_0 = [1\ 1\ 1\ 0\ 0\ 0\ 0]\). Finally, the optimal code corresponding to (63) in a matrix format is

$$\begin{aligned} \mathcal {C}_{\textrm{opt}} = \left[ \begin{array}{lllllll} 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 1 & 1 \\ 0 & 1 & 0 & 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 & 1 & 0 \end{array}\right] . \end{aligned}$$
(64)
Fig. 2
figure 2

Quantum circuit corresponding to \(E''(\textbf{x}) = \textbf{x}^{\textrm{T}} \textbf{Q} \textbf{x} + 576 - y_i\)

Full size image
$$\begin{aligned} \textbf{Q}= \left[ \begin{array}{llllllllllllllllllllll} -176 & 64 & 64 & 64 & 64 & 33 & 64 & 33 & 33 & 33 & 33 & 32 & 64 & 33 & 33 & 33 & 33 & 32 & 64 & 64 & 33 & 33 \\ 0 & -176 & 64 & 64 & 33 & 64 & 33 & 64 & 33 & 33 & 32 & 33 & 33 & 64 & 33 & 33 & 32 & 33 & 64 & 33 & 64 & 33 \\ 0 & 0 & -176 & 33 & 64 & 64 & 33 & 33 & 64 & 32 & 33 & 33 & 33 & 33 & 64 & 32 & 33 & 33 & 33 & 64 & 64 & 33 \\ 0 & 0 & 0 & -176 & 64 & 64 & 33 & 33 & 32 & 64 & 33 & 33 & 33 & 33 & 32 & 64 & 33 & 33 & 64 & 33 & 33 & 64 \\ 0 & 0 & 0 & 0 & -176 & 64 & 33 & 32 & 33 & 33 & 64 & 33 & 33 & 32 & 33 & 33 & 64 & 33 & 33 & 64 & 33 & 64 \\ 0 & 0 & 0 & 0 & 0 & -176 & 32 & 33 & 33 & 33 & 33 & 64 & 32 & 33 & 33 & 33 & 33 & 64 & 33 & 33 & 64 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 64 & 64 & 33 & 64 & 33 & 33 & 33 & 33 & 32 & 64 & 64 & 33 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 33 & 64 & 33 & 64 & 33 & 33 & 32 & 33 & 64 & 33 & 64 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 33 & 64 & 64 & 33 & 33 & 64 & 32 & 33 & 33 & 33 & 64 & 64 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 33 & 33 & 32 & 64 & 33 & 33 & 64 & 33 & 33 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 33 & 32 & 33 & 33 & 64 & 33 & 33 & 64 & 33 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 32 & 33 & 33 & 33 & 33 & 64 & 33 & 33 & 64 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 64 & 64 & 33 & 64 & 64 & 33 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 33 & 64 & 64 & 33 & 64 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 33 & 64 & 64 & 33 & 64 & 64 & 33 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 64 & 33 & 33 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 33 & 64 & 33 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 33 & 33 & 64 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 & 64 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -176 \end{array}\right] \nonumber \\ \end{aligned}$$
(65)

4 Performance comparisons

In this section, we compare the proposed GAS with the conventional GAS in terms of query complexity and demonstrate that it converges to the optimal solutions faster.

Evaluation metrics In the following, we use two metrics: (1) the total number of classical iterations required to obtain the optimal solution, which is equal to the number of measurements required, and (2) query complexity, which is equal to the total number of Grover rotations \(\sum _i L_i\). Both metrics are the same as those used in [42]. In [50], the former is called ‘query complexity in the classical domain’ and the latter is called ‘query complexity in the quantum domain.’ Typically, the query complexity in the quantum domain is considered as an important metric for evaluating quantum algorithms, corresponding to the time complexity in classical algorithms.

Simulation parameters All the parameters are the same as those used in Sect. 3.5, where we had \((n, w, d, M) = (7, 3, 4, 7)\). The objective function values were averaged over \(10^6\) trials. The termination condition of GAS can be defined by, for example, the number of iterations, time, or threshold value. Since the simulation was performed as a benchmark in this paper, we defined the termination condition that the threshold value \(y_i\) reaches the minimum objective function value. We assume the realization of FTQC.

Considered schemes The classic exhaustive search is considered as a reference scheme since the concepts of computational complexity in classical computing and query complexity in quantum computing are completely different metrics. The original search space size is calculated as \(\left( {\begin{array}{c}\left( {\begin{array}{c}n\\ w\end{array}}\right) \\ M\end{array}}\right) \). Using Theorem 3 and \(\textbf{P}'(n, w)\), the search space size can be readily reduced to \(\left( {\begin{array}{c}q_1\\ M-1\end{array}}\right) \). Here, \(q_1\), defined in (36), is always less than \(\left( {\begin{array}{c}n\\ w\end{array}}\right) \). The first codeword is fixed to \(\textbf{p}_0\), and the remaining \(M-1\) codewords are selected from \(q_1\) candidates. The conventional GAS (Algorithm 2) is considered as a performance baseline. The proposed GAS (Algorithm 3) relies on the derived bounds, which contribute to increasing the speedup. Both classical exhaustive search and conventional GAS used the proposed objective function \(E'(\textbf{x})\) of (33), while the proposed GAS used the objective function \(E''(\textbf{x})\) of (47).

For example, in the \((n, w, d, M) = (7, 3, 4, 7)\) case, the original search space size is

$$\begin{aligned} \left( {\begin{array}{c}\left( {\begin{array}{c}n\\ w\end{array}}\right) \\ M\end{array}}\right) = \left( {\begin{array}{c}\left( {\begin{array}{c}7\\ 3\end{array}}\right) \\ 7\end{array}}\right) = \left( {\begin{array}{c}35\\ 7\end{array}}\right) = 6724520. \end{aligned}$$
(66)

Here, we use Theorem 3 and \(\textbf{P}'(n, w)\). Since we have \(q_1 = 22\) from (36), the search space size is significantly reduced to

$$\begin{aligned} \left( {\begin{array}{c}q_1\\ M-1\end{array}}\right) = \left( {\begin{array}{c}22\\ 7-1\end{array}}\right) = 74613. \end{aligned}$$
(67)

The proposed formulation for GAS yields a search space size of

$$\begin{aligned} 2^{q_1} = 2^{22} = 4194304, \end{aligned}$$
(68)

and the expected minimum query complexity is

$$\begin{aligned} \sqrt{\frac{2^{q_1}}{{\underline{t}}}} = \sqrt{\frac{4194304}{6}} \sim 836. \end{aligned}$$
(69)

Figure 3 shows the expected query complexity \(k/P_k(\underline{t})\) in (56), over the interval obtained from Theorem 5. As shown in Fig. 3, when \(\overline{k}_{opt} \sim 656\), the expected query complexity is minimized, which is utilized in the proposed GAS. Note that all the following results, including conventional GAS, proposed GAS, and the classical exhaustive search, are simulated over the reduced search space of (67).

Fig. 3
figure 3

Plot of the expected query complexity \(k/P_k(\underline{t})\) in (56), over the interval obtained from Theorem 5

Full size image

First, Fig. 4 shows the relationship between each evaluation metric and the average of the objective function values, where the total numbers of classical iterations and Grover rotations were evaluated. As shown in Fig. 4, the objective function values for the conventional and proposed GAS differed significantly. By suppressing the objective function value, the number of qubits was reduced from \(q_2' = 22\) to \(q_2'' = 15\). Additionally, the proposed GAS exhibited the fastest convergence performance among the considered schemes in both metrics.

A question arises here ‘Why does the proposed GAS converge so fast?’ One major reason is the quadratic speedup of GAS with the aid of quantum superposition, entanglement, and amplitude amplification. Another reason is that the proposed GAS starts with the strict initial threshold, \(y_0 = \overline{E}'_{\textrm{min}} + 1\) of (45), and the number of Grover rotations \(L_i\) is chosen from a limited range, \([0, \bar{k}_{\textrm{opt}})\). These proposals lead to fewer states of interest and higher probability of success, resulting in faster convergence.

In Fig. 4, we simply averaged \(10^6\) curves for each scheme. After converging to the minimum objective function value, which was 15, the value of each curve remained the same. That is, in Fig. 4, the point where each scheme reached the minimum was the maximum of \(10^6\) simulations. For example, the classical exhaustive search reached the minimum objective function value at the query complexity of approximately \(5.0 \cdot 10^4\), and this complexity was the maximum of all the trials. The query complexity required to reach the minimum is an important metric and is a stochastic variable. Its statistical property is investigated in Fig. 5.

In Fig. 5, the cumulative distribution function (CDF) of the query complexity required to reach the optimum was investigated. In Fig. 4, the difference between the conventional GAS and the proposed GAS appears to be small, but when compared with CDF in Fig. 5, the difference is clear. This result indicates that on average the proposed GAS reaches the optimal solution faster than the conventional one in both the total numbers of classical iterations and Grover rotations. That is, our approach is capable of reducing the constant overhead involved in the quadratic speedup of the conventional GAS.

Fig. 4
figure 4

Relationship between each metric and the average of the objective function values, where we had \((n, w, d, M) = (7, 3, 4, 7)\)

Full size image
Fig. 5
figure 5

CDF of evaluation metrics required to reach the optimal solution, where we had \((n, w, d, M) = (7, 3, 4, 7)\)

Full size image

5 Conclusions

Although much work has been done on binary constant weight codes, its parameter-dependent combinatorial explosion is still a problem. In this paper, we proposed a quantum search algorithm for binary constant weight codes. The search problem was formulated as a polynomial binary optimization problem, and bounds were derived for the minimum of the objective function value and the exact number of solutions. We proposed a modified algorithm for GAS using the derived bounds and showed that it can further accelerate the quantum algorithm while reducing the number of qubits. The overall query complexity remains exponential because the problem is considered NP-complete. Similar approaches that utilize problem-specific bounds, such as those in this study, have the potential for wide use in applications of quantum search algorithms to other problems, thereby improving the feasibility of quantum speedup.