Finding a Hadamard Matrix by Simulated Quantum Annealing

Hard problems have recently become an important issue in computing. Various methods, including a heuristic approach that is inspired by physical phenomena, are being explored. In this paper, we propose the use of simulated quantum annealing (SQA) to find a Hadamard matrix, which is itself a hard problem. We reformulate the problem as an energy minimization of spin vectors connected by a complete graph. The computation is conducted based on a path-integral Monte-Carlo (PIMC) SQA of the spin vector system, with an applied transverse magnetic field whose strength is decreased over time. In the numerical experiments, the proposed method is employed to find low-order Hadamard matrices, including the ones that cannot be constructed trivially by the Sylvester method. The scaling property of the method and the measurement of residual energy after a sufficiently large number of iterations show that SQA outperforms simulated annealing (SA) in solving this hard problem.


Background
Finding a solution to a hard problem is a challenging task in computing. Such a problem is characterized by its complexity, as it grows beyond the polynomial against the size of the input. A class of particularly important ones are NP (non-deterministic polynomial) problems, in which verifying a solution can be conducted in polynomial time, whereas finding the solution is of exponential order. Examples of such problems are, among others, the TSP (traveling salesman problem), SAT (Boolean satisfiability), graph coloring, graph isomorphism, and subset sums.
An interesting approach to the hard problems is a method inspired by physical phenomena, such as classical annealing (CA) or quantum annealing (QA). Both CA and QA are physical processes that obtain an ordered (physical) system from an unordered one, which can be done either thermally (as is the case in CA) or quantum-mechanically (as is the case in QA). To simulate the physical processes on a (classical/non-quantum) computer, numerical methods, such as MC (Monte Carlo) for CA and PIMC (path-integral Monte Carlo) for QA, have been developed. The algorithm or computational method inspired by classical/thermal annealing is called simulated annealing (SA), whereas the one based on quantum annealing is called simulated quantum annealing (SQA). Both of these methods make use of the methods in numerical CA or numerical QA. They encode the problem into a Hamiltonian of a spin system [1] and then evolve the system from a high energy state down to the ground state. The annealing process enables the system to avoid local minima trapping and therefore is capable of achieving a global optimum, which represents the best solution of the problem. The main difference between SA and SQA is in the evolution of the systems; whereas SA uses classical/thermal annealing, SQA employs quantum mechanism.
In SA [2][3][4], one starts the system in total randomness with regard to a high temperature state. The temperature is then lowered and the system is evolved, which causes the energy to decrease so that the system becomes increasingly ordered. To avoid local-optima trapping, a particular updating rule, such as the Metropolis [2], is applied. The rule allows the system to (sometimes) move to a higher energy state. Upon completion of the algorithm, the system achieves the ground state, at which point a solution is found.
In [5], Kadowaki and Nishimori introduced quantum fluctuations to replace the thermal fluctuations in SA to accelerate the convergence. They applied the method on an Ising model, where a transverse field plays the role of temperature in classical SA, enabling the system to achieve the ground state with greater probability. Santoro et al. [6] compared classical and quantum Monte Carlo annealing protocols on a two-dimensional Ising model. They found that the quantum Monte Carlo annealing is superior to classical annealing. In [7], Boixo et al. show experimental results on a 108 qubit D-Wave One, which is a kind of hardware implementation of QA. A strong correlation between D-Wave and SQA, compared to the device with classical annealing, was found, which indicates that the D-Wave performs quantum annealing. This result raised the important issue of whether QA actually outperforms SA [8]. Rønnow et al. [9] showed how quantum speedup should be defined and measured. In an experiment with random spin glass instances on 503 qubits of D-Wave Two, they did not find any evidence of such speedup.
Regardless of these issues, different results have been achieved via SQA. Isakov [10] performed quantum Monte Carlo (QMC) simulations and found that the QMC tunneling rate displayed scaled according to system size. He also found quadratic speedup in QMC simulations when, instead of periodic conditions, open boundary conditions were employed. In [11], Mazzola et al. demonstrated that QMC simulations can recover the scaling of ground-state tunneling rates, which validates QA in terms of solving combinatorial problems.
Some classes of hard problems, including ones with exponential or combinatorial complexity, have been a subject of interest in SQA research. Martonak et al. [12] introduced an application of SQA to solve the TSP problem. They found that a PIMC algorithm was more efficient than SA in terms of finding an approximately minimal tour in a given graph. SQA has also been used to successfully address other hard problems related to graphs, such as graph coloring [13] and graph isomorphism [14].
In this paper, we propose SQA as a mean to find a Hadamard matrix (H-matrix). Previously, in [15], we successfully employed SA to perform a similar task, in which low-order H-matrices were found. Compared to existing H-matrix construction methods, an SA-based method is more general in terms of its capability of finding (or constructing probabilistically) an m = 4k order H-matrix, without any restriction on the property of the order m, whereas the Sylvester method requires m = 2 n , where k and n are positive integers. This paper extends this classical SA method to its quantum version, where PIMC based on Suzuki-Trotter formulation [16,17] is employed to simulate the quantum process.

Finding A Hadamard Matrix
A Hadamard matrix, or H-matrix, is an orthogonal binary {±1} matrix of size 4k × 4k, where k is a positive integer. This matrix was discovered by J. J. Sylvester [18] in 1867 and then studied more extensively by J. Hadamard [19] during his investigation of the maximal determinant problem. The orthogonal property makes the H-matrix popular in applied areas, such as information coding and signal transform. In the 1960s and 1970s, Hadamard code was used in space exploration for information transmission [20,21]. In a recent technological case, CDMA (code-division multiple access), which is widely used in cellular mobile phone systems, employs Walsh-Hadamard signals to reduce interference between its users [22,23].
One of the most important issues in the theory of H-matrix is its existence. Any 2 l order H-matrix with l a positive integer can be constructed using Sylvester's method. Furthermore, if there is an m order H-matrix, m = 4k can be shown for a positive integer k. On the other hand, no one yet knows if there is always a 4k order H-matrix [20,21]. The latter case is formulated as the Hadamard matrix conjecture. Up to this writing, the smallest unknown 4k order H-matrix is 668.
Various reconstruction methods have been proposed [24][25][26][27][28][29]. Nevertheless, these methods force the order m to follow a particular rule. In [15], a general m = 4k order algorithm employing SA is proposed. The method works on a special H-matrix called a seminormalized Hadamard (SH) matrix, in which the first column is a 4k order unity vector v 0 = (1, · · · , 1) T , and the rest are 4k order SH vectors v i ∈ V.
A brute-force method needs to verify all N B of the 4k order binary matrix to find an H-matrix, where N B (4k) = 2 16k 2 [15]. Let all matrices constructed where v 0 is the first column and a combination of v i ∈ V constitutes the remaining (4k − 1) columns be called quasi-SH (QSH) matrices. Since there are unique QSH-matrices. Although the number has been greatly reduced compared to N B , exhaustive checking still requires a great amount of computational resources. The SA method proposed in [15] is capable of finding a few low-order SH matrices in a more reasonable time.
Following the convention in our previous paper [15], the role of the spin, i.e., its ±1 eigenvalues, is replaced by SH spin vectors v i ∈ V. To find a 4k order SH-matrix, one needs (4k − 1) fully connected SH spin vectors, which initially are set randomly. With a defined energy E( Q), the SH spin vectors are randomly changed in accordance with conditions whereby a transition into another SH spin vector is allowed but a transition into a non-SH-spin-vector is forbidden.

Simulated Quantum Annealing
The Hamiltonian of an Ising system with spin configuration {σ k }, where k ∈ K = {1, 2, · · · , i, j, · · · } is the set of the lattice's indices, can be expressed aŝ where J ij is a coupling constant/strength between a spin at site i with a spin at site j, h j is the magnetic strength at site j, and {σ z i ,σ x i } are Pauli's matrices at site i. In SQA, quantum fluctuation is elaborated by introducing a transverse magnetic field Γ. The Hamiltonian of the system takes the following form [5]:Ĥ In Equation (2), the transverse field is changed (reduced) over time, i.e., Γ ≡ Γ(t). On the right hand side of the equation, the first two terms corresponds to potential energyĤ pot , while the third one is the Hamiltonian introduced by the transverse field, which is related to kinetic energyĤ kin ; i.e, we can defineĤ In general,Ĥ pot andĤ kin do not commute, so [Ĥ pot ,Ĥ kin ] = 0. Denoting the Hamiltonian of the potential as a function of spin configurationsĤ pot ≡Ĥ {σ z i } , we can also express Equation (2) in a more general form as follows:Ĥ To simulate a quantum system described by Equation (5) using the classical method, we have to formulate PIMC by introducing imaginary time. It can be then approximated by the Suzuki-Trotter transform by adding one dimension in the imaginary time direction, which, for (P × N) degrees of freedom, takes the following form [13,30]: where N is the number of spins in the lattice, P is the number of Trotter's replicas, S i = ±1 are the eigenvalues of the spin matrices, and is the nearest-neighbor coupling of the transverse magnetic field [30].

SQA Formulation of the SH Spin Vector
Similar to the previous paper [15], we employ a seminormalized Hadamard spin vector, abbreviated here as an SH spin vector, instead of an ordinary spin. In a 4k order SH spin vector, for a given positive integer k, 2k spins are −1 and another 2k spins are +1. Therefore, an SH spin vector transition is allowed only if these balance numbers are conserved; otherwise, such a transition is forbidden. We also treat the SH spin vector as a single entity, even though it consists of 4k spins, and is denoted as v i ∈ V, where V is the set of all 4k-order SH vectors. We formulate the energy of a particular configuration of spin vectors { v i } as follows: where v i · v j denotes the inner product of the vector v i with v j . Figure 1 shows an Ising system with four SH spin vectors with an additional Trotter's dimension. In the lower part of Figure 1a, each circle represents a binary spin, whereas the solid line represents the connection among the spins. Interacting spin i with binary variable S i and spin j with binary variable S j contributes the term J ij S i S j to the Hamiltonian. For a 4k order case, every 4k non-connected spins are grouped into one SH vector v i , which is illustrated as a dashed line. To simplify the diagram, each SH vector is represented by a filled circle; thus, we obtain the upper part of Figure 1a, which is called a slice or a replica. In the PIMC, the slice is replicated P-times, and these slices are arranged as layers in imaginary time. Each neighboring SH vector in a replica, i.e., v i,p with v i,p−1 and v i,p with v i,p+1 , interacts. The extension (in imaginary time) is illustrated in Figure 1b. The Hamiltonian in Equation (6) becomes a Hamiltonian of an SH vector spin system H QV that can be rewritten as follows: where J Γ ≡ J Γ (t) and H pot { v i,p } represent complete-graph connections among the SH spin vectors, similar to Equation (8), which is given by The evolution of H QV in Equation (9)  We will now formulate the SQA method for finding the H-matrix into an algorithm, which is displayed as pseudo-code in Algorithm 1. It takes the matrix order, the number of replicas, the initial temperature, the initial value of Γ, and the amount of iterations and sub-iterations as inputs. This algorithm yields either an SH-matrix or a QSH-matrix that has more orthogonal column vectors than the initial one. The algorithm starts with a random initialization of replicas with QSH-matrices, which are (4k − 1) sets of SH vectors, and then calculates its initial energy. Following the schedule of a linear transverse field, a trial transition is performed for each replica. The acceptance and rejection of the transition is based on the Metropolis criterion. The iteration will be stopped when either the number of maximum iterations is reached or an SH-matrix is found.

Finding a 12-Order SH-Matrix Using SQA
We have performed numerical experiments to find low-order H-matrices. Here we present results for the H-matrix of order 12 for detailed analysis, since it is the lowest-order H-matrix that cannot be constructed by the Sylvester method. Initially, all of the slices (replica) were filled with randomly generated v i ∈ V. Note that there are two nested iterations in Algorithm 1. The first one is an iteration of all replicas with the maximum number set to k · M × M, where M = 12 is the H-order. The second one is an iteration of flipping within a slice of a replica, whose number is c · M, c can be any small number.
The energy evolution during the iteration is shown in Figure 2. The figure shows curves of replica energy E rep , mean potential energy Ep mean , and minimum potential energy Ep min . The replica energy is defined similarly to Equation (9), i.e., E rep ≡ H QV , whereas the potential energy is given by Equation (10) E pot ≡ H pot . The mean and minimum values have been taken across the replicas. Based on the figure, both Ep mean and E rep fluctuate over time, but they tend to decrease. The minimum energy of a lattice in the replica Ep min also tends to decrease. When Ep min = 0, the H-matrix is found. Figure 2. Energy evolution during the SQA algorithm runs to find an SH-matrix of order 12. Four curves are drawn in the graph, which are the mean potential energy Ep mean , the minimum potential energy Ep min , the replica energy E rep , and the deviation standard of the potential energy Ep std . When Ep min equals zero, the iteration is stopped since an SH-matrix has been found. The Ep std curve indicates high variation in the configuration of replicas at the initial stage, which is then reduced in later stages.
The degree of orthogonality of the matrix Q is displayed by the indicator matrix D ≡ Q T Q. Figure 3 shows the initial QSH-matrix and its related indicator matrix. We also show the initial and final indicators for the first and last slices of the replica in Figure 4. It is expected that all of the QSH-matrices become more orthogonal, indicated by a lower number of zeros in off-diagonal entries. The last figure showing the last slice of the replica condition after the iterations are completed clearly show this case. The found H-matrix is shown in the left part of Figure 5, with its corresponding indicator shown on the right, which is a diagonal matrix.

The Number of The Replicas and Convergence Issue
In theory, the number of Trotter's replicas P should be as large as possible. However, in practice, we should also consider the convergence issue when a running time restriction (iteration number) is given. As explained in [13,30], replicas provide diversity of solutions; a greater P selects the best solution with minimum energy. On the other hand, the replicas are not merely running Monte Carlo on several replicas; the interactions between replicas J Γ (t) also define their behavior, i.e., a large value of Γ at the initial stage implies a low value of J Γ , which loosens the connections, and the interactions then become independent. A low Γ value at the end of an iteration implies a high J Γ value, which tightens the replica connections such that they become similar. To measure these variations, we used a simple deviation standard of energy across the replicas. Figure 6 shows the curves of variation of the energy evolution for P = 5, 10, 15, and 20 in finding a 12-order H-matrix. Figure 6. The effect of the replica number P in the algorithm: although ideally a large P is desired, it also needs to be adjusted to the problem. Variation in replica energy (in terms of deviation standards of the energy across the replicas) when searching for an H-matrix is shown. The numbers of replicas P = 10, 15, 20 yield large variations up to the end of the iteration, whereas P = 5 yields a better result with steady values at the end. In all of these cases, for the construction of a 12-order H-matrix, the total maximum iteration is set to 20,000, consisting of a global iteration count of 20,000 for each P.
Since initially the replicas were set randomly, they will have almost identical energy, so variation in the energy will be very low. In later iterations, the value will increase as a new configuration is explored, and this will be followed by a decrease, which indicates that the replicas have become homogeneous. This cycle of increasing-decreasing energy should be observed if P is chosen properly with respect to the dimension of the problem (H-order) and a sufficient number of iterations. When P is too small, the system will perform akin to classical SA, whereas a P that is too large will cause the system to fail. The figure shows that, for a given number of maximum iterations 20,000, the number of replicas P = 5 is the most suitable; anything higher is too high. This also shows that frequent updates on a limited number of replicas, compared to less frequent updates on a larger number of replicas, better achieve convergence.

Performance Comparison: SQA vs. SA
To compare performances, in the first experiment, we measured the residual error of both algorithms. Since the ground state is achieved when the matrix becomes orthogonal, in which case Equation (10) will equal zero, the residual error will be defined as the minimum H pot over all of the replicas, i.e., = min (H pot ). We have chosen the order of the H-matrix to be sufficiently large so that we will still have a residual error at the end of the execution of the algorithm, i.e., so that the H-matrix is not found. We considered order M = 28 to be sufficient for this purpose, where we actually have 28 3 = 21,952 spins. We also chose a Trotter slice of P = 5 and plotted the curve for iterations 50 up to 5,000,000.
Following [30], the annealing schedule was linear; i.e., the temperature T was reduced linearly in SA, and was the transverse magnetic strength Γ. Even though T is reduced linearly, the threshold probability P thresh will change exponentially. By using the function the threshold will start a bit higher than 0.5, which asymptotically approaches 1.0 at the end of iteration time t. Figure 7 shows the curve of T(t), P thresh (t), Γ(t), and J Γ (t). The experiments were repeated 10 times for each case. The averages of residual errors for each iteration numbers are plotted in Figure 8 for both SA and SQA.
The figure shows that, although initially the residual error of SQA is larger than SA, the slope is steeper. With a higher number of iterations, which in this case is around 100,000, SQA is superior. Considering that SQA shows the least amount of error among the replica slices, it seems that variation in the replica is an ideal solution. In SA, once a solution is selected, the change in spin configuration will be less significant by the time the system reaches a lower energy state. Therefore, in terms of finding an H-matrix, SQA is superior to SA. In the second experiment, both SA and SQA were applied to matrices with an increasing size (order). Figure 9 shows a graph of computational gain, which is defined as the ratio of the number of SA iterations to the number of SQA iterations needed to achieve 50 percent of the residual energy of the initial mean energy of all replicas. The horizontal axis shows the order of the H-matrix, from 4 to 20, whereas the vertical axis shows the computational gain. The gain grows with the order of the H-matrix, which shows that speedup increases with problem size. Based on this curve, we observe that SQA outperforms SA for the Hadamard search problem. Figure 9. Curve of computational gain, which is the ratio of the number of SA iterations to the number of SQA iterations needed for the algorithm to achieve 50 percent of its initial residual error. The horizontal axis represent the problem size, which is the order of the H-matrix. The figure shows that the gain grows non-linearly with problem size, indicating that SQA outperforms SA.

Conclusions
We here propose a new method of finding an H-matrix based on SQA. We have formulated the method into an algorithm, which has been implemented, tested, and analyzed. Low-order H-matrices, including one of order 12 that cannot be constructed via the Sylvester method, were found. We have also discussed the advantages of the method over classical SA. Measurements of the residual error and the relative running time on an increasing order of H-matrices indicate that SQA is superior to SA in solving the Hadamard search problem.