Quantum computing formulation of some classical Hadamard matrix searching methods and its implementation on a quantum computer

Finding a Hadamard matrix (H-matrix) among all possible binary matrices of corresponding order is a hard problem that can be solved by a quantum computer. Due to the limitation on the number of qubits and connections in current quantum processors, only low order H-matrix search of orders 2 and 4 were implementable by previous method. In this paper, we show that by adopting classical searching techniques of the H-matrices, we can formulate new quantum computing methods for finding higher order ones. We present some results of finding H-matrices of order up to more than one hundred and a prototypical experiment of the classical-quantum resource balancing method that yields a 92-order H-matrix previously found by Jet Propulsion Laboratory researchers in 1961 using a mainframe computer. Since the exactness of the solutions can be verified by an orthogonality test performed in polynomial time; which is untypical for optimization of hard problems, the proposed method can potentially be used for demonstrating practical quantum supremacy in the near future.


Quantum computing formulation of some classical Hadamard matrix searching methods and its implementation on a quantum computer Andriyan Bayu Suksmono 1* & Yuichiro Minato 2
Finding a Hadamard matrix (H-matrix) among all possible binary matrices of corresponding order is a hard problem that can be solved by a quantum computer. Due to the limitation on the number of qubits and connections in current quantum processors, only low order H-matrix search of orders 2 and 4 were implementable by previous method. In this paper, we show that by adopting classical searching techniques of the H-matrices, we can formulate new quantum computing methods for finding higher order ones. We present some results of finding H-matrices of order up to more than one hundred and a prototypical experiment of the classical-quantum resource balancing method that yields a 92-order H-matrix previously found by Jet Propulsion Laboratory researchers in 1961 using a mainframe computer. Since the exactness of the solutions can be verified by an orthogonality test performed in polynomial time; which is untypical for optimization of hard problems, the proposed method can potentially be used for demonstrating practical quantum supremacy in the near future. Background. A Hadamard matrix (H-matrix) is a binary orthogonal matrix with {−1, +1} elements whose any distinct pair of its columns (or rows) are orthogonal to each other. Such a matrix only exists when it is square and the length of its column (row) is equal to 1, 2 or a multiple of four; i.e., for an M × M dimension H-matrix, then M = 1, 2 or M = 4k for a positive integer k. The reversed statement that for any positive integer k there is a H-matrix is also believed to be true, although neither a mathematical proof nor disproof yet exists. This is a long standing problem of the Hadamard Matrix Conjecture.
The H-matrix has been a subject of scientific and practical interests. First discovered and described by Sylverster in 1867 1 , it is further studied by Hadamard concerning its relationship with the determinant problem 2 . The orthogonal property and binaryness of its elements make it widely used in information processing and digital communications. The CDMA (Code Division Multiple Access) system employs Hadamard-Walsh code to reduce interference among their users, so that the capacity of the communication system is not badly deteriorated by the increasing number of its users 3,4 . The H-matrix was also used by Mariner 9 space-craft as its ECC (Error Correcting Code) for sending images of Mars to a receiving station located on Earth, thanks to its capability for long error correction 4,5 .
Some particular kinds of H-matrices can be found (constructed) easily, while others need huge computational resource to do. An H-matrix of size M × M is also called an M-order H-matrix. When M follows a particular pattern of M = 2 n , where n is a positive integer, the matrix can be easily constructed by the Sylvester's method of tensor product. Hadamard 2 constructed the H-matrices of order 12 and 20, whose orders do not follow the 2 n pattern. It indicates that other orders than prescribed by the Sylvester's method do exist. Paley showed the construction of H-matrix of order M = 4k where k ≡ 1 mod 4 and k ≡ 3 mod 4 , which are known as the Paley Type I and Type II H-matrices, respectively 6 . In the formulation, he employed the method of quadratic residues in a Galois field GF(q), where q is a power of an odd prime number.
Various kinds of construction methods have also been proposed. A cocyclic technique, which is based on a group development over a finite group G modified by the action of a cocycle defined on G × G , has been OPEN 1 The School of Electrical Engineering and Informatics, Institut Teknologi Bandung, Bandung, Indonesia. 2  www.nature.com/scientificreports/ introduced by De Launey and Horadam 7,8 . The Hadamard matrices can be generated by this scheme when it is applied to binary matrices. A general introduction on the cocyclic methods are described by Horadam 9 and recent progress are presented; among others by Alvarez et al. 10,11 . In developing a quantum computing based H-matrix searching method, we found that a simple and straight forward method will be a good starting point. Our methods described in this paper have been based on earlier techniques proposed by Williamson 12 , Baumert-Hall 13 , and Turyn 14 , which is suitable for this purpose. These three methods involve searching of particular binary sequences as an essential stage. In this paper, we will refer these methods to as classical H-matrix searching methods. Although at a glance it looks simple, finding a H-matrix is actually a challenging task. To find a H-matrix of order 92, in 1961 three JPL (Jet Propulsion Laboratory) researchers employed a state-of -the-art computer at that time, i.e. the IBM/7090 Mainframe 15 . For matrix order under 1000, the most recent unknown H-matrix successfully found is the one with order 428, which was discovered in 2005 by using computer search of particular binary sequences 16 . The method described in the paper is of particular interest because the next unknown H-matrices, such as the one with order 668, possibly can be found by using the same method. The main reason they have not been found at this time is because of the huge computational resource needed to find such matrices, which grows exponentially by the order of the matrix.
Finding a H-matrix of order M among all of O(2 M 2 ) binary matrices, which we refer to as H-SEARCH, is a hard problem. We have proposed to find such a matrix by using a quantum computer considering its capability in solving hard problems 17 . Theoretically, a quantum computer will need O(M 2 ) qubits in superposition to solve such a problem. However, in the existing quantum annealing processor, we need O(M 3 ) due to extra ancillary qubits required to translate k-body terms into 2-body Ising Hamiltonian model. In this paper, we show that by adopting the classical searching methods, we can reduce the required computing resource, which for a quantum annealing processor implementing the Ising model, will become O(M 2 ) . We describe how to formulate the corresponding Hamiltonians related to the classical methods and show some results of order up to more than one hundred. We also describe how to further develop this technique to find higher order matrices, by managing the classical and quantum computing resources. In such a classical-quantum hybridized algorithm, the complexity of the classical part still grows exponentially, but the quantum part grows polynomially. We shows that this algorithm extends the capability of a pure quantum method with limited number of qubits, so that a few higher order of H-matrices can be found, compared to the pure quantum method that cannot be implemented on present days quantum computer.
Usually, solving an optimization problem by annealing or heuristic methods yields only an approximate solution, i.e., we can not sure that it is actually the optimal point, unless all of possible solutions are enumerated. However, enumeration of all possible solutions of a hard problem is an extremely laborious task. In contrast, the correctness of a solution in H-SEARCH can be verified easily in polynomial time; i.e., by evaluating the orthogonality of the found matrix (solution). If we consider the solution as a certificate, H-SEARCH behaves like an NP-complete problem because finding the solution is hard, but checking its correctness is easy. In this particular point of view, H-SEARCH is an interesting hard problem worth to consider in addressing practical quantum supremacy.
A brief on quantum computing and finding H-matrices using quantum computers. Quantum computers are expected to have computational capability beyond their classical counterparts; a feature which is well known as quantum speedup 18 or even quantum supremacy 19 . An important progress regarding this issue is the achievement of the Google researchers in 2019, who claimed that their Sycamore quantum processor needs only about 200 s to do a particular computational task; which is sampling random quantum circuits in this case, where a classical supercomputer would take about 10,000 years to perform 20 . In the next step, a capability of solving a real-life problems, where classical computers cannot do in a reasonable time, is desired. Creative thinking of building algorithms that can demonstrate such practical supremacy are needed.
The working principle of QAM computers are based on quantum annealing (QA) 22,23 , which is a quantum analog to the classical (thermal) annealing (CA). Whereas the CA works by gradually decreasing temperature with sometimes allowing the system to jump over higher energy, the QA seeks for the solution by quantum tunneling through the energy barrier. Energy landscape of the H-SEARCH problem's Hamiltonian are degenerates; i.e., there are many equivalent binary matrices that have identical energy. Illustration of the potential energy landscapes (PEL) for 2-order and 4-order binary matrices are given in Fig. 1. Considering the PELs, quantum annealing approach is suitable to find the solution. We expect that the speed-up comes from the process of finding the minimum energy by quantum tunneling. Further analysis on how quantum computing can speed up a search algorithm is described by Farhi and Gutmann 24 . A comprehensive review on quantum annealing and analog quantum computation has been given by Das and Chakrabarti 25 .
In general, existing quantum computers can be categorized into the universal quantum gate (QGM-Quantum Gate Machine) and quantum annealer (QAM-Quantum Annealing Machine). Regardless some issues related to noise and other non-ideal conditions, both of these types of quantum processors have been built and are accessible by public users through the Internet. The implementation scheme of the proposed methods for both of these kinds of quantum computers are illustrated in Fig. 2. The direct method; which works for QAM that has been described in our previous paper 17 , will be used as a reference. Three main proposed quantum computing methods are derived from non-quantum computing/classical H-matrix construction/searching methods, which we will referred to as the Williamson, Baumert-Hall, and Turyn methods.
The QAM processor, such as the D-Wave, only accepts problems in the form of a 2-body Hamiltonian that generally can be expressed by www.nature.com/scientificreports/ which is a Hamiltonian of an Ising system, where J ij is a coupling constant or interaction strength between a spin at site i with a spin at site j, h j is magnetic strength at site j, and {σ α i } are Pauli's matrices of directions α = {x, y, z} at site-i. The processor performs quantum annealing by introducing a transverse field given by which is evolved over time according to the following equation where t ∈ [0, τ ] denotes time 22,26 . The problem to solve should be encoded in Ĥ pot , which is represented by the Ising's coefficient J ij and h i for each of the problem. Some optimization problems have been solved by the quantum annealing methods; among others are: graph isomorphism 27 , wireless network optimization 28 , nurse scheduling problem 29 , hand written digit recognition 30 , computational biology 31 , and hydrologic inverse analysis 32 .
In a QAM, the formulation of the H-SEARCH is started by calculation of its energy function E(s) as a function of binary variables s ∈ {−1, +1} . For conciseness, we will represent the value of s by its signs {−, +} . In general, E(s) might contain high order k-body interaction terms so that we will denote it by E k (s) , whereas the Ising model allows only up to 2-body terms in E 2 (s) . To obtain the 2-body expression, and eventually a 2-body quantum Hamiltonian Ĥ 2 σ , a sequence of transforms given by the following construction diagram should be conducted 17 , where q ∈ {0, 1} is a Boolean variable. Actually both of s and q are binary variables, but with different values. For now on, we will refer s ∈ {−1, +1} as spin variable and q ∈ {0, 1} as Boolean variable.
In the previous paper 17 , implementation of an M-order H-matrix on a QAM needs M 2 number of logical (binary) variables and additional M 2 × (M − 1)/2 ancillary variables (ancillas) so that the overall complexity is O(M 3 ) . In this paper, by adopting classical H-matrix construction/searching methods, we can reduce the required number of variables significantly into O(M 2 ) which enables the search of higher order H-matrices than before. In the followings, we will address three quantum H-SEARCH methods, which are derived from the classical methods of Williamson, Baumert-Hall, and Turyn. For each of these methods, we derive their corresponding Hamiltonians based on some criteria that are specifics for each of the cases. Low order cases can be calculated by hand, while higher order ones should be calculated by a computer through symbolic computing due to the  21 . In (a), the 2-order binary matrices are generated and their energies are plotted against their indices. The indices are converted to the matrices after a binary to spin variable transform. As an example, for the second matrix with index 1, then the process is as follows: In the QGM quantum computing, we can employ QAOA (Quantum Approximate Optimization Algorithm) 33 , which is well-suited for solving an optimization problem on NISQ (Noisy Intermediate-Scale Quantum) processors. In principle, the general k-body Hamiltonian can directly be implemented on a QGM. Therefore, the required number of physical qubits will be about the same as the number of logical qubits. However, since the implementation needs direct connection to the actual machine, which is not available for us at this time, we will not address it in the current paper.

Results
Williamson based quantum computing method. The Williamson's method builds a matrix W of size 4k × 4k from four sub-matrices A, B, C, D each of size k × k 4,12,34 . Any pair of these sub-matrices are commutative. The orthogonality property of W will be satisfied when where V = I k is a k × k identity matrix. Then, the problem becomes choosing the elements of s i ∈ {−1, +1} in the sub-matrices that makes the orthogonality condition in Eq. (5) is satisfied. Further simplification and efficiency of the number of variables can be achieved when we choose the sub-matrices which are symmetric and circular.
By imposing the orthogonality conditions, the commutativity among the sub-matrices, and the non-negativity of the energy, we arrive to the following s-dependent energy function where v i,j denotes the element at row i and column j of the matrix V that consists of products of spin/binary variables s i given by Eq. (5) and δ i,j is the Kronecker delta function. The orthogonality requirement of W will be satisfied when E k (s) = 0 , which is the lowest value of the energy function of Eq. (6). For k = 3 and a particular www.nature.com/scientificreports/ value of Boolean reduction factor δ (note that it was written as δ ij in 17 ), by expanding this equation and then following the construction diagram in Eq. (4), we will arrive to the following 2-body Hamiltonian which can be encoded into a quantum annealing processor.
In the experiment, we extract the Ising coefficients {J ij , h i } then submit them to the D-Wave. We observe that the magnitude of the coefficients in the Hamiltonian's terms are quite large, however they will be normalized by the D-Wave system. Additionally, the constant term, such as 162, 720 in Ĥ 2 (σ z ) of Eq. (7), will also be removed. Consequently, instead of zero, the minimum of the energy will be a negative value. We have set the number of reads to 10,000 and obtain some solutions at minimum energy values. For k = 3 , which corresponds to H-matrix of order 12, the required number of logical qubits was 8 which translates into 36 physical qubits. We obtained the minimum energy at −45.988 . The experimental results are displayed in Fig. 3a, where the bottom part shows the found H-matrix H on the left side and its indicator matrix D ≡ H T H on the right side, whereas the top parts show energy distribution of the solutions. Higher orders matrices, up to order 36 that needs 49 physical qubits to implement, have also been found successfully using the D-Wave. They are listed in the Supplementary Information section.
Baumert-Hall based quantum computing method. The Baumert-Hall method works in a similar manner as the Williamson's by first finding the A, B, C, D block matrices, except that the construction of the H-matrix is given by a 12 × 12 structure of block matrix 13,34 , which yields a 12k × 12k matrix for particular values of positive integers k.
Experiments on finding Baumert-Hall matrices using D-Wave quantum processor indicates that the capability of the method is limited by the available number of qubits, the number of couplers, and the capability of the embedding tool 35 . We have successfully found a few of Hadamard matrices up to order 108 using this method. For the 108-order case; which corresponds to k = 9 , by following the construction diagram with particular value of the Boolean reduction factor δ , we will obtain a 2-body Hamiltonian given by, After extracting the Ising parameters and submitting to the D-Wave, we obtain the solutions containing correct values of s i for building the H-matrices. Figure 3b shows a 108 order H-matrix, which was found by the Baumert-Hall based method and its corresponding energy statistics as output of the quantum computer. Other Baumert-Hall matrices found by this method, i.e. 36, 60 and 84, are listed in the Supplementary Information section.
The Turyn based quantum computing Method. In this method, first we have to find a set of 4-sequences {X, Y , Z, W} that has particular properties, then use them to construct a H-matrix based on Goethals-Seidel method 14,16 . We derive the energy function from the requirement of a valid TT-sequences given by, where N X (r), N Y (r), N Z (r), N W (r) are non-periodic auto-correlation functions of the sequences {X, Y , Z, W} calculated at lag-r, respectively. Since the value given by the left-hand side of Eq. (9) can be negative, whereas the annealing is performed to achieve a minimum value, we modify it into a non-negative energy function which are squared sum of the auto-correlation function at each lag r ≥ 1 as follows, We have inserted a k subscript to indicate that the energy may includes k-body interaction terms. The searching problem becomes finding a TT-sequence that satisfy this condition. We will represent the elements of {X, Y , Z, W} as spin variables s i as before. As an example we will calculate the Hamiltonian for k = 4 . By considering normalized sequence for efficiently use variables 16 , we obtain the following expressions for TT (4) Then, the energy in Eq. (10) which after following construction diagram given by Eq. (4) with a particular value of Boolean reduction factor δ , yields the following 2-body Hamiltonian  TT-sequence needs 2 4n−1 steps, which is an exponentially increasing task. For finding higher order H-matrices, we can explore the properties of the TT-sequence to reduce the number of binary sequence to enumerate 16,36 .
In this method, instead of finding all {X, Y , Z, W} at once, it will be more computationally realistic to start with filling some part of them, then subsequently imposing conditions and properties of the TT-sequence to limit the number of the sequences to check.
Partially filled sequences {X * , Y * , Z * , W * } with m-elements on the left part and another m-elements on the right one, are given as follows The requirement of non-periodic auto-correlation sum for these sequences is now become We will refer all {X * , Y * , Z * , W * } sequences satisfying condition given by Eq. (13) as solution prototypes. Then the energy function becomes  www.nature.com/scientificreports/ TT-sequence of up to 40 can be found using classical computers, whereas higher order ones need more powerful computers which is impossible to be implemented at the moment. This is one of the main reasons that H-matrix of order 668 has not been found nor declared non-exists yet, assuming that such a matrix can be constructed by the Turyn's method.
On the other hand, we can use the solution prototypes to reduce the number of required qubits when a quantum computer is involved in the searching process. For clarity, in the followings we illustrate this method by a simple case which is implementable on a current quantum processor. We will consider a (4, 4, 4, 3) solution prototype to find a (8,8,8,7) TT-sequences by using quantum computing; therefore, it is a kind of finding higher order sequence by extending the lower one. The extended TT-sequences can be expressed by with known x 0 , . . . , x 3 , y 0 , . . . , y 3 , . . . , . . . , w 0 , w 2 and unknown s 0 , s 1 , . . . , s 15 .
To find the unknown values represented by s i , we calculate the energy of the Turyn's based method as before. Among all possible {X * , Y * , Z * , W * } prototypes and the replacement of the unknowns with binary variables, we choose the following solution prototype as an example Note that in the real case, we might have to check all of the solution prototypes.
Further calculation by applying the construction diagram with a particular value of Boolean reduction factor δ yields the following 2-body Hamiltonian, We then encode the Hamiltonian into the D-Wave. Calculation shows that we need 108 physical qubits to implement, but embedding into the Chimera graph with the D-Wave provided embedding tools indicates that more qubits are required, which in this case is 860. After quantum annealing, we get among others, the following solution Among 10,000 of obtained results, we identified two correct solutions. One of the solution that has been constructed to a H-matrix; its corresponding indicator matrix, and solution statistics are displayed in Fig. 3d. This TT(8)-sequences yields a 92-order Hadamard matrix, which in 1961 was also found by JPL researchers in a search using IBM/7090 mainframe computer 15 .

Discussions
Difficulties in finding a H-matrix by classical computing methods, due to the exponential grows of the complexity, can be overcome by quantum-computing based search such as by directly represents each elements of the matrix into binary variables, which is then translated into qubits 17 . However, the availability of quantum computing resource limits the implementation to only finding low order H-matrices. We have shown in the previous section that classical H-matrix searching methods can be adopted to efficiently use available quantum computing resource to solve larger problems, i.e., finding higher order H-matrices than previously achieved by the direct method 17 .
The data displayed in the top part of Table 1 shows required resource and results in the Williamson and Baumert-Hall based methods for each order of the H-matrix. Since both of them share the same A, B, C, D block matrices, we put them side-by-side on the table. We observe in the table that the number of required logical qubits grows linearly by O(M) with respect to the order of the searched matrix, whereas the number of physical qubits grows quadratically as O(M 2 ) , which is caused by the ancillary qubits required to translate k-body into 2-body Hamiltonians. In the implementation, the physical qubits and their connections should be mapped to the topology of qubits's connections in the quantum annealing processor; which is the Chimera graph in the DW-2000Q. We have used (default) embedding tool provided by D-Wave 35 and the number of embedding qubits displayed in the table are taken from the output of the software. This mapping process, which is also called minor www.nature.com/scientificreports/ embedding, further increases the number of required qubits. In the following discussions, the number of required qubits after the embedding process will be labeled as the embedding qubits. The Williamson and Baumert-Hall adopted methods can be implemented to all of matrix order as long as the embedding process is successful, which is up to 36 for the Williamson and up to 108 for the Baumert-Hall. We observe from the output of embedding tool that the highest order needs 1, 492 qubits to implement, which is more than 6 times of the required physical qubits. Observing that the trends of the embedding-to-physical qubits ratio; denoted by E/P ratio in the table, increases with the H-matrix order, by taking a moderate estimate of 7 times, the 300 physical qubits for the order of 132 matrix (in the Baumert-Hall based method) requires 2, 100 qubits to be implemented; which is more than currently available qubits in the DW-2000Q. We also observe from the experiment results, especially those displayed in the last column of Table 1, that the number of correct solutions among 10,000 number of reads tends to decrease with the increasing order of the matrix; i.e., it is about 4% at the beginning then decreased to about 0.2% at order 108 for the Baumert-Hall. A possible explanation to this phenomena is that when the order of the matrix is increased, the magnitude of the coefficients in the Hamiltonian are also increased so that the difference between the largest to the smallest value becomes very large. Since they are normalized when fed into the D-Wave, with limited resolution to 1/256, the D-wave will set lower coefficients to zero. Accordingly, some of the terms will be discarded and the solutions become degenerate. It makes the Table 1. Resource required in the proposed quantum computing methods. The Willamson and Baumert-Hall based method: to find an M order H-matrix, the required logical qubits grows by O(M) and the physical qubits by O(M 2 ) . Embedding the connections implied by the Hamiltonian on existing Chimera graph further increases the required qubits, which ultimately limit the capability of the method. Decreasing percentage of correct solutions, knowing that only 10,000 reads in a single run is allowed, indicates that repeated experiments will be needed to find higher order matrices. The Turyn-based quantum computing method: although the number of required physical qubits also grows by O(M 2 ) , the jump among the order is high so that the next one after order 68, i.e. 92 and beyond,cannot successfully be embedded in DW-2000Q. We also cannot conclude the success rate for given 10,000 number of reads due to lack of data, although we may suspect that it will also decrease as in the Williamsons and Baumert-Hall adopted methods. Extended Turyn based method: the number of logical qubits is determined by the number of additional k in the extension k , not by the final qubits. We show the result of extending order k = 4 into k = 8 and only one of successful solution prototypes in the  Since the number of reads in one run is limited by the D-Wave system to 10,000, several repeated runs on the quantum processor should be done to find higher order H-matrices. Figure 5 plots the probability of success against the order of Baumert-Hall matrices; it shows that in general higher order matrices are more difficult than the lower ones to find by the method. This also means that, for finding higher order H-matrices; assuming that the processor has a number of sufficient qubits to implement, what we have to do is by repeating the experiments many times.
Middle part of Table 1 shows required number of qubits and performance of the Turyn based quantum computing method. An H-matrix of order 44 and 68 have successfully been found, but higher order matrices have not. We observe that the E/P ratio grows faster than the similar case in the Williamson and Baumert-Hall based methods; i.e, it is now about 11 times at the order of 68. Assuming this factor stay the same, higher order matrices of 92, which needs 199 physical qubits, might require about 2, 189 embedding qubits. This is more than the currently available number of qubits in the DW-2000Q quantum processor, and therefore the search of order 92 H-matrix did not successful. We have proposed a solution for the limitation of quantum computer resource by the extended Turyn based method described previously.
Bottom part of Table 1 shows the required resource and performance of extended Turyn based method. For extending k = 4 into k = 8 , we need 108 physical qubits; which is then increased into 860 embedding qubits. An important feature of this method is that, as long as the number of additional/extension k = 4 is kept the same, the required qubits to solve extended problem will also stay the same, regardless the targeted order. However, this advantage should be paid by increasing number of solution prototypes, implying that more classical computing resources is needed and the frequency usage of the quantum processor will be increased. We expect to have an optimal point where the combination of the classical and quantum resources delivers the best solution and achieves highest order of the searched H-matrix.
At present, some of H-matrices of order under 1000; such as 668, 716, and 892, have not been found by any methods yet due to huge computational resource required to perform the computation by existing classical methods. When using the Turyn-based quantum computing method, even after extension, H-SEARCH for such orders still cannot be implemented. As an illustration, with a 12 pre-filled {X * , Y * , Z * , W * } , the required logical qubits for the 668 case will be 4 · (56 − 12) = 176 which becomes 176 · 175/2 = 15, 400 physical qubits. Assuming the similar embedding performance as before at a factor of 8, the required number of qubits is 123, 200 which is beyond the capability of current quantum annealing processors. Figure 6 shows the progress of available qubits in D-Wave quantum annealers [? ] and the decrease number of required qubits to implement H-matrix search of order 668 by solving the {X * , Y * , Z * , W * } problem. The points in the graph shows actual number of qubits achieved in every year since 2011. We can see that the number of qubits doubled every two years; therefore, by using regression we get a linear line in a semi-logarithmic plot as shown by a dotted red curve. The middle dashed green horizontal line indicates the number of required qubits when no additional embedding qubits are required, which means that an ideal complete graph connection among the qubits is available. The top blue dashed dotted line indicates the number of required qubits with embedding factor of 8. Assuming that the connections among qubits are also improved substantially every year, we can expect the H-SEARCH of order 668 can be implemented between the year 2022 to 2029. Additionally, recent achievement of 64 qubits volume 43 and the 1000 qubits milestone 44 of QGM processor, the H-SEARCH implementation through QAOA is also very promising to explore. where H 2 = + + + − , i.e., it is a kind of plugging-in smaller H-matrices into a particular structure to obtain a larger H-matrix. Similarly, the Williamson method also builds a higher-order matrix from smaller ones, except that the smaller matrices are not necessarily orthogonal. In general, we can express the Williamson type H-matrices W by 4,12,34 where A, B, C, D are block matrices, whose any pair of them are commutative, i.e., where δ i,j is the Kronecker delta function. The orthogonality requirement of W will be satisfied when E k (s) = 0 , which is the lowest value of the energy function in Eq. (23). In the k = 3 case, the energy function E k (s) can be expanded into For implementing an energy function to a QAM processor; such as in the D-Wave, the k-body energy function E k (s) should be transformed into a 2-body energy function E 2 (s) using the steps given by the construction diagram in Eq. (4). In the process, we should choose a Boolean reduction factor δ to transform the k-body into 2-body function, that should be larger than the maximum value E max of the energy function 45 . By taking E max = 26, 976 , which is the maximum value of E k (s) by assuming all of s i = +1 , then setting δ = 2E max , we obtain the following result This 2-body energy function gives the potential Hamiltonian Ĥ pot (σ ) ≡Ĥ 2 (σ z ) as follows, Complete expressions for the equations can be found in Supplementary Information section.

(22)
24) E k (s) = 6(4 + 2(s 0 s 1 + s 2 s 3 + s 4 s 5 + s 6 s 7 )) 2 (25) E 2 (s) = 13, 728s 0 + 13, 728s 1 + · · · + 13, 488s 0 s 1 + · · · + 192s 10 s 11 + 162, 720 Experiments on finding Baumert-Hall matrices using D-Wave quantum processor indicates that the capability of the method is limited by the available number of qubits and the capability of the embedding tool 35 . We have successfully find the Hadamard matrix up to order 108 using this method. For the 108-order case, initial energy function E k (s) to find this matrix is given by the following whose corresponding k-body Hamiltonian is given by Then the 2-body Hamiltonian realized on the quantum annealing processor will be given by, Complete expressions for the equations can be found in Supplementary Information section. Derivation of the Turyn based method. In this method, first we find a set of 4-sequences {X, Y , Z, W} that has particular properties, then use them to construct a H-matrix based on Goethals-Seidel method 14,16 . We translate the requirements into energy functions which then programmed into a quantum processor. In essence, the workflows of the Turyn based method are as follows 1. Find an (n, n, n, n − 1) Turyn-Type (TT) sequence {X, Y , Z, W}. 2. Construct base sequences {A, B, C, D} 3. Construct T-sequences {T 1 , T 2 , T 3 , T 4 } 4. Construct seed sequences {A 1 , A 2 , A 3 , A 4 } 5. Construct block symmetric circular matrices {X A1 , X A2 , X A3 , X A4 } 6. Construct Hadamard matrix, which is given by where R is a back-diagonal identity matrix of size k × k as follows We derive the energy function from the requirement of a valid TT-sequences given by, E k (s) = 432s 0 s 2 + · · · + 720s 18 s 19 + · · · + 432s 16 s 17 s 18 s 19 + 5760 0 0 · · · 0 1 0 0 · · · 1 0 · · · · · · · · · · · · · · · 0 1 · · · 0 0 1 0 · · · 0 0      www.nature.com/scientificreports/ where N X (r), N Y (r), N Z (r), N W (r) are non-periodic auto-correlation functions of the sequences {X, Y , Z, W} calculated at lag-r, respectively. The non-periodic auto-correlation function of a sequence V = [v 0 , v 1 , . . . , v n−1 ] T is given by, for r = 0, 1, . . . , n − 1 and N V (r) = 0 for r ≥ n . Since the value given by the left-hand side of Eq. (33) can be negative, whereas the annealing is performed to achieve a minimum value, we adopt a non-negative energy function which are sum of squared value of the auto-correlation function at each lag r ≥ 1 as follows, To efficiently use available qubits in the quantum processor, it is important to reduce the number of variables encoded to the qubits as few as possible. We can achieve this by further employing the property of the TT-sequences. In this case, we can normalize the TT-sequences 16 to obtain X T = (x 0 , x 1 , . . . x n−1 ) , Y T = (y 0 , y 1 , . . . y n−1 ) T , Z T = (z 0 , z 1 , . . . z n−1 ) , and W T = (w 0 , w 1 , . . . w n−1 ) , which have the following properties For clarity, in the followings we present an example of the Hamiltonian formulation for the lowest order of k = 4 case. The first step as described previously is to find a TT(4) -sequences {X, Y , Z, W} . By representing the elements of the sequences as binary (spin) variables s i ∈ {−1, +1} , and applying the properties of a normalized sequence explained previously, a TT(4) will be as follows, To determine the energy function, we have to calculate non-periodic auto-correlation functions N X , N Y , N Z , N W given by Eq. (34). Since s 2 i = 1 , we get the following results after simplifications Therefore, the energy E ≡ E k (s) in Eq. (35), whose terms may contain products of k variables, is now given by In the following steps, as described by the construction diagram in Eq. (4), the energy function should be transformed into a 2-body interacting Ising Hamiltonian. Therefore, we have to change the s-dependent energy function into q-dependent energy function E k (q) . After simplification, this transform yields the following The conversion into 2-body energy function requires a Boolean reduction factor δ set to be larger than the maximum value of the energy function E max (k) . Assuming it is at least an absolute sum of the E k (q) coefficients as before, we have E max = 908 . By taking twice of this maximum value, we obtain δ = 1, 816 , which transforms Eq. (39) into E k (s) = 2s 1 + 2s 2 + 2s 4 + 4s 0 s 3 + 4s 1 s 2 − 4s 0 s 4 + 2s 1 s 3 + 2s 1 s 4 + 2s 2 s 3 + 2s 2 s 4 + 4s 0 s 1 s 2 + 2s 1 s 2 s 3 + 4s 0 s 3 s 4 + 2s 1 s 3 s 4 + 2s 2 s 3 s 4 + 2s 1 s 2 s 3 s 4 + 242 E k (q) = −16q 0 − 40q 1 − 40q 2 − 40q 3 − 24q 4 + 16q 0 q 1 + 16q 0 q 2 + 32q 0 q 3 + 48q 1 q 2 + 32q 1 q 3 + 24q 1 q 4 + 32q 2 q 3 + 24q 2 q 4 + 40q 3 q 4 − 32q 0 q 1 q 2 − 32q 1 q 2 q 3 − 32q 0 q 3 q 4 − 16q 1 q 2 q 4 − 32q 1 q 3 q 4 − 32q 2 q 3 q 4 + 32q 1 q 2 q 3 q 4 + 276  Further reduction of the needed qubits is achieved through the usage of proposed methods described in this paper, such as the Turyn based method. As explained in the section Methods, subsection C. Derivation of the Turyn based Method, the (Turyn) Hadamard matrix can be constructed from an (n, n, n, n − 1) Turyn Type/ TT-sequence. For a given (n, n, n, n − 1) TT-sequence, we can construct a 4(4n − 1) order Hadamard matrix; i.e, to find an 4(4n − 1) order H matrix, we only need to find a (4n − 1) length sequence. Therefore, in the Turynbased method, the required logical qubits to find the M × M Hadamard matrix is O(M). The quadratic energy function given by Eq. (35) implies that there will be up to 4-body terms in the Hamiltonian. Again, when using D-Wave that can only accommodate up-to 2-body terms, additional ancillary qubits are needed. Accordingly, the final number of the required logical qubits will be O(M 2 ).

Data and code availability
All of codes and data will be provided upon direct request to the authors. Some parts of the codes can be found in https:// github. com/ suksm ono.