Dual-Matrix Domain-Wall: A Novel Technique for Generating Permutations by QUBO and Ising Models with Quadratic Sizes

The Ising model is defined by an objective function using a quadratic formula of qubit variables. The problem of an Ising model aims to determine the qubit values of the variables that minimize the objective function, and many optimization problems can be reduced to this problem. In this paper, we focus on optimization problems related to permutations, where the goal is to find the optimal permutation out of the $n!$ possible permutations of $n$ elements. To represent these problems as Ising models, a commonly employed approach is to use a kernel that utilizes one-hot encoding to find any one of the $n!$ permutations as the optimal solution. However, this kernel contains a large number of quadratic terms and high absolute coefficient values. The main contribution of this paper is the introduction of a novel permutation encoding technique called dual-matrix domain-wall, which significantly reduces the number of quadratic terms and the maximum absolute coefficient values in the kernel. Surprisingly, our dual-matrix domain-wall encoding reduces the quadratic term count and maximum absolute coefficient values from $n^3-n^2$ and $2n-4$ to $6n^2-12n+4$ and $2$, respectively. We also demonstrate the applicability of our encoding technique to partial permutations and Quadratic Unconstrained Binary Optimization (QUBO) models. Furthermore, we discuss a family of permutation problems that can be efficiently implemented using Ising/QUBO models with our dual-matrix domain-wall encoding.


Introduction
A Binary Quadratic Model (BQM) [1] is defined by an objective function that includes a quadratic formula with multiple variables.The problem associated with a BQM is to find the values of these variables that minimize the resulting value of the quadratic formula.A BQM is referred to as a Quadratic Unconstrained Binary Optimization (QUBO) [2] model when the variables are restricted to binary values, i.e., they can only take bit (or binary) values in {0, 1}.On the other hand, if the variables can only take qubit (or spin) values in {−1, +1}, the model is called an Ising model.It is worth noting that QUBO and Ising models can be converted into each other interchangeably [1,3].
Since optimization problems such as the traveling salesman problem; scheduling problems; and various graph problems, including max cut, maximum independent set, and graph isomorphism, can be transformed into QUBO/Ising models [4], there has been significant research dedicated to finding efficient algorithms, hardware, and systems to solve them.However, solving optimization problems for QUBO/Ising models is known to be NP-hard.
This means that unless P = NP, it is not possible to design a polynomial time algorithm using classical computers with digital circuit devices of polynomial size.In the quest for solutions, researchers have explored the potential of ideal quantum annealers based on quantum mechanics, which could potentially find optimal solutions for large Ising models within a reasonable time frame [5].Unfortunately, the current quantum annealers available are not yet powerful enough to tackle such problems effectively.The number of qubits is limited, and the presence of undesirable flux noise significantly reduces the probability of finding optimal solutions [6].Hence, as an alternative to an ideal quantum annealer, BQM solvers on various non-quantum computing platforms, such as ASICs [7,8], FPGAs [9][10][11][12], GPUs [13][14][15][16], and optimal fibers [17,18], have been proposed.Further, D-Wave Systems released a hybrid BQM solver [19] that uses both a classical computer and a quantum annealer to find solutions for large BQMs with up to 1,000,000-node complete graphs.
In this paper, our primary focus is on the size and resolution of QUBO/Ising models, which are represented by the quadratic term count and the maximum absolute value of the coefficients, respectively.These metrics serve as important indicators for assessing the models.While it is not always the case, QUBO/Ising models with smaller sizes and resolutions are generally considered more desirable.This is particularly true when considering the limitations of quantum annealers, which have restricted size and resolution capabilities.Therefore, it becomes crucial to ensure that the Ising models embedded in quantum annealers possess small sizes and low resolutions.Additionally, even when these models are processed by classical digital computers, the required memory size to store QUBO/Ising models is proportional to the quadratic term count, and a higher resolution requires a larger word size of memory.With this perspective in mind, this paper places significant emphasis on the size and resolution of QUBO/Ising models.We provide precise evaluations of these metrics, offering valuable insights into their characteristics.
QUBO/Ising models that are converted from permutation-based combinatorial optimization problems for n elements should have a kernel capable of generating any one of the n! possible permutations as the optimal solution.To achieve this, a common approach is to use the one-hot encoding of permutations, which involves a bit/qubit matrix of size n × n.In this encoding, each row i (0 ≤ i ≤ n − 1) of the matrix represents the i-th number as a one-hot vector, where exactly one element is set to 1/+1, and its position corresponds to the number it represents.The reader should refer to Figure 1, which illustrates a 4 × 4 matrix representing the permutation [1,3,2,0].In order to generate any one of these matrices as the optimal solution, QUBO/Ising models require kernels that have optimal solutions if and only if each row and each column contains exactly one 1/+1.For instance, QUBO models utilizing this kernel have been proposed for addressing Hamiltonian cycle/path problems and the graph isomorphism problem, as demonstrated in [4].Similarly, models for the traveling salesman problem (TSP) and the quadratic assignment problem (QAP) were introduced in [20,21].However, these kernels typically involve n 3 − n 2 quadratic terms.Additionally, the kernel of the Ising model features a non-constant large coefficient of 2n − 4. Recently, a permutation generation technique using domain-wall encoding has been introduced [22].This technique employs a matrix of size n × (n − 1), where each row i (0 ≤ i ≤ n − 1) stores the i-th number as a domain-wall vector [23][24][25].By utilizing this technique, the number of quadratic terms is reduced to 1  2 n 3 − 3 2 n, which is half the number of quadratic terms in QUBO/Ising kernels obtained through conventional one-hot encoding.This technique still involves a cubic number of quadratic terms.Moreover, QUBO/Ising models using this domain-wall encoding require coefficients with large absolute values, namely 2n − 3/n − 1, respectively.
The main contribution of this paper is to introduce a novel permutation encoding technique called the dual-matrix domain wall, which uses two matrices A and B to store dual permutations.This new technique can significantly reduce the number of quadratic terms and the maximum absolute coefficient values in the resulting Ising kernel.The QUBO/Ising kernel obtained by this technique has only 6n 2 − 12n + 4 quadratic terms.Also, the maximum absolute value of the coefficients of the Ising kernel is just 2. We also discuss partial-permutation generation by QUBO/Ising kernels, which is a sequence of m numbers selected from n numbers 0, 1, . .., n − 1 without repetition.It is not difficult to apply one-hot encoding to generate a partial permutation [26,27].For partialpermutation generation, one uses an m × n matrix such that each i-row (0 ≤ i ≤ n − 1) stores the i-th selected number as a one-hot vector.QUBO/Ising kernels that use one-hot encoding to generate any such matrices need 1  2 m 2 n + 1 2 mn 2 − mn quadratic terms.We show that we can apply the dual-matrix domain-wall technique for partial-permutation generation, and the resulting QUBO/Ising models have only 6mn − 6m − 6n + 4 quadratic terms.
Moreover, we introduce a generic problem called the particle placement problem (PPP), which aims to optimize the placement of m particles in n locations (m ≤ n) with no collision.In other words, the problem is to find an optimal permutation of m numbers selected from n numbers.In this study, we demonstrate that the PPP can be efficiently reduced to a QUBO/Ising model using permutation-generating kernels.Additionally, we establish that several permutation-based combinatorial optimization problems, such as the quadratic assignment problem (QAP), the traveling salesman problem (TSP), the sub-graph isomorphism problem, and the maximum weight matching problem, can also be reduced to the PPP.Thus, these problems can be converted to equivalent QUBO/Ising models using permutation-generating kernels.
Finally, we conducted an evaluation of the cells required to embed Ising kernels on a D-Wave quantum annealer, specifically the Advantage 4.1.The objective was twofold: first, to assess the feasibility of utilizing the dual-matrix domain-wall technique on presently accessible quantum annealers, and second, to examine the influence of the number of quadratic terms in Ising models on the cell requirements.
This paper is organized as follows.In Section 2, we begin by providing a formal definition of QUBO and Ising models.We also establish the relationship between quantum annealers and Ising models.Additionally, we explore different encoding techniques such as one-hot, zero-one-hot, and domain-wall encoding for representing numbers in QUBO and Ising models.Section 3 reviews the conventional one-hot encoding technique used to represent permutations of n numbers.Furthermore, we introduce the all-different domain-wall encoding technique, which effectively reduces the number of quadratic terms by half [22].However, both of these techniques require a cubic number of quadratic terms.In Section 4, we present a novel encoding technique called dual-matrix domain-wall encoding.This technique enables the generation of permutations using a QUBO/Ising kernel with only a quadratic number of quadratic terms.Section 5 generalizes the concept of permutations to partial permutations, which represent a partial permutation of m numbers selected from a set of n numbers without repetition.We demonstrate how QUBO/Ising kernels can generate partial permutations using both the conventional one-hot encoding technique and our dual-matrix domain-wall encoding technique.In Section 6, we extend the dual-matrix domain-wall technique by incorporating a matrix called a one-hot matrix that stores a partial permutation as a one-hot encoding.We introduce the particle placement problem (PPP) in Section 7, highlighting its ability to reduce many permutation-based combinatorial optimization problems.We evaluate the number of quadratic terms in the resulting QUBO/Ising models obtained through reduction via the PPP for each problem.Section 8 presents an analysis of the number of cells required to embed Ising kernels for generating permutations in a D-Wave quantum annealer.Finally, Section 9 concludes our work.

Preliminaries
The main objective of this section is to introduce the concepts of QUBO and Ising models, as well as the fundamental aspects of their design.To begin with, we provide a formal definition of Ising and QUBO models.Subsequently, we establish the connections between Ising models and quantum annealers, offering valuable insights to the reader.Lastly, we showcase the applications of QUBO and Ising models in generating one-hot, zero-one-hot, and domain-wall vectors.

QUBO and Ising Models
A Binary Quadratic Model (BQM) [1] is defined as an objective function with a quadratic formula involving multiple variables.The goal of a BQM is to determine the variable values that minimize the resulting value of the quadratic formula.
A model is referred to as a Quadratic Unconstrained Binary Optimization (QUBO) model [2] if the variables can take bit (or binary) values of 0 or 1.Specifically, let X = (x i ) (0 ≤ i ≤ n − 1) represent an n-bit vector of binary variables.A QUBO model with X can be defined using an upper triangular weight matrix W = (W i,j ) (0 ≤ i ≤ j ≤ n − 1).The objective of the QUBO problem is to find the binary values of X that minimize the energy E(X), which is defined by the following quadratic formula: Here, C is a constant called the offset.Note that Equations ( 1) and ( 2) are equivalent because x 2 i = x i always holds.While most papers use the energy E(X) in Equation ( 2) with no offset C, this paper adopts the energy E(X) defined by Equation (1) to distinguish linear and quadratic terms with coefficients W i,i (0 ≤ i ≤ n − 1) and W i,j (0 ≤ i < j ≤ n − 1), respectively.
A BQM is called an Ising model if variables take qubit (or spin) values of −1 or +1.An Ising model with an n-qubit vector S = (s i ) (0 ≤ i ≤ n − 1) is defined by quadratic term coefficients J = (J i,j ) (0 ≤ i < j ≤ n − 1) called interactions and linear term coefficients h = (h i ) (0 ≤ i ≤ n − 1) called biases.The Ising problem aims to find the qubit values of S that minimize the Hamiltonian H(S), which is defined by the following quadratic formula: We make the assumption that all coefficients in the linear and quadratic terms of QUBO/Ising models are integers, unless otherwise specified.Using integers with large absolute values can lead to discrete values below the effective resolution, so it is important to keep the maximum absolute value of all coefficients as small as possible.To achieve this, we design QUBO/Ising models that have no common factor in all coefficients.This allows us to reduce the coefficients by dividing the common factor without affecting the optimal solutions.It is worth mentioning that the offset C does not influence the optimal solution and can be disregarded.Consequently, it can take a non-integer value.
We define the number of non-zero elements in W i,i /h i (0 ≤ i ≤ n − 1) as the linear term count and that in W i,j /J i,j (0 ≤ i < j ≤ n − 1) as the quadratic term count of a QUBO/Ising model.Generally, having smaller term counts is advantageous as it helps reduce the hardware resource usage of quantum computers when solving problems of QUBO/Ising models.In particular, the quadratic term count has a significant impact on the hardware resource usage of quantum annealers when solving Ising problems.By minimizing the quadratic term counts, we can effectively reduce the amount of resources needed for implementing and solving QUBO/Ising models on quantum hardware.This optimization is essential for improving the efficiency and performance of quantum computing systems.
It is easy to see that QUBO and Ising models can be equivalently converted to each other [3].Both QUBO and Ising models, excluding the offset C, can be represented as weighted undirected graphs with n nodes 0, 1, . . ., n − 1 such that W i,i /h i is the weight of node i and W i,j /J i,j is the weight of edge (i, j) of the graph.Figure 2 shows a graph corresponding to a QUBO model with the the following energy: The figure also illustrates the equivalent Ising model with the following Hamiltonian: We omit terms with zero coefficients in the formulas and edges with zero weights in the graphs.As QUBO and Ising models are represented as graphs, we use graph theory terminologies such as the degree of a node (i.e., the number of edges connecting to a node) and the diameter (i.e., the largest shortest path over all pairs of nodes).The QUBO and Ising models have optimal solutions X = [1, 0, 0, 1, 1, 0] with energy E(X) = −7 and S = [+1, −1, −1, +1, +1, −1] with Hamiltonian H(S) = −32, respectively.These models are equivalent, because 4E(X) = H(S) + 4 always holds for all X and S satisfying s i = 2x i − 1 for all i (0 ≤ i ≤ n − 1).As many optimization problems can be reduced to QUBO/Ising models [4], there has been a significant effort by researchers to find effective algorithms, hardware, and systems to solve them.Digital computers and devices can directly operate 0/1 bits, and users/developers can handle them more easily than −1/+1 qubits.As a result, QUBO models are more frequently used than Ising models, and QUBO solvers have been developed by many researchers.On the other hand, quantum annealers based on quantum mechanics can directly operate −1/+1 qubits in Ising models.Therefore, for solving QUBO problems on quantum annealers, they are converted to equivalent Ising models.By this preprocessing step, quantum annealers can be used as QUBO solvers.However, to maximize the performance of QUBO solvers, we may design Ising models without using QUBO-Ising conversion.
When a QUBO/Ising model for solving a specific permutation-based combinatorial optimization problem is designed, this involves creating a sub-model that generates any one of all possible permutations as the optimal solution.This sub-model is referred to as the QUBO/Ising kernel.The main focus of this paper is the design of QUBO/Ising kernels for generating permutations.For instance, we will present QUBO/Ising kernels for generating an n × n-matrix that stores one-hot vectors representing a permutation, as shown in Figure 1.

Quantum Annealers and Ising Models
D-Wave Systems developed a quantum annealer called D-Wave 2000Q [28], which was based on quantum mechanics.It served as a solver for Ising models using a 2048-node Chimera graph.Subsequently, D-Wave Systems released a more advanced quantum annealer known as D-Wave Advantage [29], which was capable of handling Ising models with a larger 5760-node Pegasus graph [30].Specifically, the D-Wave Advantage quantum annealer comprises 5760 cells interconnected according to the topology of a 5760-node Pegasus graph.The cell biases and interaction strengths are programmable, allowing for the acquisition of optimal or approximate solutions to the corresponding Ising models through quantum annealing.
While a quantum annealer is designed to solve Ising models with a specific graph topology, it is possible to solve Ising models with different topologies through a process known as minor embedding.Minor embedding involves mapping the problem graph (with a different topology) onto the physical graph of the quantum annealer, effectively embedding the problem into the hardware.This allows the quantum annealer to solve Ising models that may not directly match its native graph topology.Minor embedding is a technique used to leverage the capabilities of quantum annealers for a broader range of Ising models.We explain the idea of minor embedding using Figure 3. Suppose that we need to solve an Ising problem with the six-node graph in Figure 2 on a quantum annealer with an eight-cell grid topology, as shown in Figure 3.We arrange node 0 (or qubit s 0 ) in Figure 2 to two cells 0 (or qubit s 0 ) and 0 ′ (or qubit s 0 ′ ) in Figure 3.We add a quadratic term −Ps 0 s 0 ′ so that optimal solutions satisfy s 0 = s 0 ′ , where P is a constant number large enough to guarantee it.By virtue of this embedding, s 0 in Figure 2 can be simulated by s 0 and s 0 ′ combined in Figure 3.The reader should have no difficulty in confirming that the Ising model in Figure 2 can be solved by the quantum annealer in Figure 3.More than two cells in the quantum annealer may be arranged to simulate a node of the Ising model.D-Wave Systems calls a set of cells arranged to form a qubit a chain, and the value of P the chain strength [31].To obtain an optimal or approximate solution of an Ising model by the D-Wave Advantage quantum annealer, one needs to find its minor embedding in a 5760-node Pegasus graph.If an Ising model is a dense graph with a large degree, the chains will be large.For example, if an Ising model is a complete graph, only 177 nodes can be embedded in the D-Wave Advantage quantum annealer [29].Hence, the size of an Ising model that can be solved by the D-Wave Advantage quantum annealer is limited.In particular, the quadratic term count of an Ising model impacts the resource usage of cells.Thus, it is quite important to minimize the quadratic term count when one designs Ising models.
The D-Wave Advantage quantum annealer [29] can take real numbers for the interaction J i,j and bias h i of Ising models.The ranges of J i,j and h i are limited to [−1.0, +1.0] and [−4.0, +4.0], respectively.Although they take any real number in these ranges, they are operated as analog values, and the effective resolution is limited.The resolution is only five to six bits, and two values with a difference less than 1  2 6 = 0.015 may not be distinguished due to the flux noise.Recall that we assume the interaction J i,j and bias h i of Ising models are integers.When we load such Ising models onto the D-Wave Advantage quantum annealer, these integer values are automatically reduced to fit in the ranges [−1.0, +1.0] and [−4.0, +4.0].For example, if J i,j takes integers in the range [−100, +100], then they will be divided by 100 to fit in the range [−1.0, +1.0], and the required resolution 0.01 is smaller than the effective resolution of the D-Wave Advantage quantum annealer.Hence, optimal solutions of the Ising model cannot be expected.Since the resolution is limited, the maximum absolute values of interactions and biases of Ising models should be as small as possible to obtain optimal solutions with a higher probability by quantum annealers.
QUBO/Ising models with larger diameters may be solvable using a heuristic algorithm in parallel.They may be split into many disjoint sub-models, and parallel divide-andconquer techniques can be applied to such models.Heuristic algorithms to improve the solutions of sub-models can be executed in parallel, and better approximate solutions can be obtained more quickly.On the other hand, smaller-diameter models cannot be split into many disjoint sub-models, and it is difficult to apply parallel heuristic techniques.Thus, larger-diameter models are preferable to solve the problems in parallel.
2.3.QUBO/Ising Models for One-Hot/Zero-One-Hot/Domain-Wall Vectors One-hot encoding has often been used to represent integers in QUBO and Ising models.A k-bit/qubit vector is a one-hot vector if it has exactly one 1/+1, and the remaining k − 1 bits/qubits take 0/−1.It represents an integer i (0 ≤ i ≤ k − 1) if and only if the i-th bit/qubit is 1/+1.For technical reasons, we introduce the zero-one-hot vector, which can take a bit/qubit vector with all 0/−1s in addition to one-hot vectors.Such vectors with all 0/−1s represent a special value φ associated with "undefined" or "N/A".Table 1 shows four-bit/qubit one-hot/zero-one-hot vectors with corresponding integers or φ.

One-Hot
Zero-One-Hot Domain-Wall Bits Qubits Bits Qubits Bits Qubits Domain-wall encoding [22][23][24][25] has also been used to represent integers.A k-bit/qubit vector is a domain-wall vector if it consists of consecutive 1/+1s following consecutive 0/−1s.It represents the integer i (0 ≤ i ≤ k) if it contains i consecutive 1s.Therefore, it can represent k + 1 integers ranging from 0 to k. Table 1 illustrates four-bit/qubit domain-wall vectors that represent integers from 0 to 4.
This sub-section presents QUBO/Ising models for one-hot/zero-one-hot/domain-wall vectors, aimed at helping the reader understand the fundamental concepts of using these models to represent numbers.These models are designed to achieve their optimal value if and only if the vectors X/S are one-hot/zero-one-hot/domain-wall vectors.
We will first design a QUBO model with k-bit variable X = (x i ) (0 ≤ i ≤ k − 1) that takes the minimum value of 0 if and only if it stores a k-bit one-hot vector.It is clear that X is a one-hot vector if and only if the sum of all bits is 1.Based on this property, we can design a QUBO model with an energy function E 1 (X) as follows: It is evident that E 1 (X) takes the minimum value of 0 if and only if X is a one-hot vector.Additionally, we can design a QUBO model that takes the minimum value if and only if it stores a k-bit zero-one-hot vector.A k-bit vector X is a zero-one-hot vector if and only if the sum of all bits is either 0 or 1.Thus, the energy E 01 (X) defined below takes the minimum value of 0 if X is a zero-one-hot vector.
QUBO models with E 01 (X) for zero-one-hot vectors have no linear terms, making them simpler than those with E 1 (X) for one-hot vectors.Since all quadratic terms in the expanded summation are 2, E 01 (X) has a coefficient of 1  2 .Our objective is to minimize the integer coefficients of both linear and quadratic terms.Consequently, we will employ such coefficients whenever possible throughout the remainder of this paper.
Similarly, we can design Ising models with H 1 (S) and H 01 (S) for a k-qubit variable S = (s i ) (0 ≤ i ≤ k − 1) storing one-hot and zero-one-hot vectors, respectively.A k-qubit variable S stores a one-hot vector if and only if Thus, the following H 1 (S) takes the minimum value of 0 if and only if S is a one-hot vector: Clearly, H 1 (S) takes the minimum value of 0 if and only if S is a one-hot vector.Additionally, S stores a zero-one-hot vector if and only if it has zero or one +1.Thus, the following Ising model with H 01 (S) takes the minimum value of 0 if and only if S is a zero-one-hot vector, which means that the sum of all qubits is −k or −(k − 2): We will now design a QUBO model for a k-bit variable X = (x i ) (0 ≤ i ≤ k − 1) that stores a domain-wall vector.For technical reasons, we assume the existence of fixed guard bits x −1 = 1 and x k = 0 for X.With these guard bits, if X stores a domain-wall vector, then x i−1 − x i ̸ = 0 (0 ≤ i ≤ k) holds for exactly one i.Otherwise, it holds for more than two i.Based on this fact, we have the following QUBO model: The energy function E d (X) takes the minimum value of 1 2 if and only if (x i−1 − x i ) 2 = 1 for exactly one i and X stores a domain-wall vector.Similarly, we can design an Ising model for a k-qubit variable S = (s i ) (0 ≤ i ≤ k − 1) that stores a domain-wall vector using the same approach.We also assume the existence of fixed guard qubits s −1 = +1 and s k = −1.
The following Ising model with H d (S) takes the minimum value of 2 if and only if S stores a domain-wall vector and (s i−1 − s i ) 2 = 4 for exactly one i: Table 2 summarizes the features of QUBO/Ising models with k-bit/qubit one-hot/zeroone-hot/domain-wall vectors.We utilized SymPy [32], a Python library for symbolic mathematics, to expand the mathematical formulas and derive the features of the QUBO/Ising models presented throughout this paper.We observe that QUBO/Ising models for one-hot vectors are fully connected and consist of 1  2 k 2 − 1 2 k quadratic terms.In contrast, models for domain-wall vectors have only k − 1 quadratic terms.Furthermore, the linear term coefficients of Ising models for one-hot and zero-one-hot vectors are k − 2 and k − 1, respectively.On the other hand, those for domain-wall vectors have only two linear terms with coefficients −1 and +1, respectively.Table 2. QUBO/Ising models for k-bit/qubit one-hot/domain-wall vectors.

Encoding Type
One-Hot Zero-One-Hot Domain-Wall

QUBO models
Quadratic formula Optimal energy 0 0 Ising models Quadratic formula

QUBO/Ising Model Kernels for Generating Permutation Involving a Cubic Number of Quadratic Terms
A permutation of n numbers can be defined by a bijection π : {0, 1, . . ., n − 1} → {0, 1, . . ., n − 1}, where the list [π(0), π(1), . . ., π(n − 1)] represents one of the n! possible permutations.This section initially describes a conventional permutation encoding method that employs one-hot vectors.This approach is widely used for solving permutation-based combinatorial optimization problems, not just in QUBO/Ising models but also in mixedinteger programming.We then explain a permutation encoding method using domain-wall vectors, which was introduced in [22] and reduces the number of quadratic terms by half.

Conventional Permutation Encoding by One-Hot Vectors
A permutation π is commonly represented by an n × n matrix of variables, where each row i (0 ≤ i ≤ n − 1) stores π(i) as a one-hot vector.This representation ensures that each row of the matrix is a one-hot vector.We will design n 2 -bit QUBO/Ising kernels with X = (x i,j )/S = (s i,j ) (0 ≤ i, j ≤ n − 1), which can generate permutations as the optimal solutions.For this purpose, we apply E 1 (X)/H 1 (S) to all rows and columns to guarantee that they are one-hot vectors as follows: We can obtain the quadratic formulas for QUBO and Ising models by expanding E nn 1 (X) and H nn 1 (S), respectively.Both formulas take the minimum value of 0 if and only if every row and every column is a one-hot vector.The features of QUBO/Ising models obtained by expanding E nn 1 (X)/H nn 1 (S) can be found in Table 3.It is important to note that both models have n 3 − n 2 quadratic terms, and the Ising model also includes linear terms with a coefficient of 2n − 4.
Table 3. QUBO/Ising kernels for generating a permutation of n numbers.

Domain-Wall Encoding Type
One-Hot All-Different Dual-Matrix Extended Ising models Quadratic formula

Permutation Encoding by All-Different Domain-Wall Encoding
We will now introduce the all-different domain-wall technique, which was presented in [22] and can be used to generate permutations with domain-wall vectors.This technique utilizes a matrix X = (x i,j ) of size n × (n − 1), where each row i (0 ≤ i ≤ n − 1) stores the domain-wall vector representing π(i) for a permutation π. Figure 1 illustrates an example with n = 4.It is important to note that a matrix X, where each row contains a domain-wall vector, represents a permutation if and only if all the domain-wall vectors are distinct.The all-different domain-wall technique leverages this property to design a QUBO kernel for generating a permutation, as follows: Here, fixed guard bits x i,−1 = 1 and x i,n−1 = 0 for all i (0 ≤ i ≤ n − 1) are used.The first summation term takes the minimum value of 1  2 n if and only if every row is a domain-wall vector.Also, the the second summation term takes the minimum value of 0 if and only if the number of 1s in each column j (0 ≤ j ≤ n − 1) is n − j − 1.If this is the case, all rows of X store distinct domain-wall vectors.Therefore, X is a permutation if and only if it is an optimal solution that satisfies E nn a (X) = 1 2 n.We can apply the same technique for designing an Ising kernel with an n × (n − 1) matrix S = (s i,j ) that generates a permutation as a domain-wall vector.If all rows in S store distinct domain-wall vectors, then the number of +1s in each column j (0 ≤ j ≤ n − 1) is (n − j − 1), and the sum of all elements in it is (n Based on this fact, we can design a desired Ising model as follows: Here, s i,−1 = +1 and s i,n−1 = −1 for all i (0 ≤ i ≤ n − 1).The first summation term takes the minimum value of 1  2 • 4n = 2n if and only if every row is a domain-wall vector.Also, the second summation term takes 0 if and only if the number of +1s in each column j (0 ≤ j ≤ n − 1) is n − j − 1.If this is the case, all rows of S store distinct domain-wall vectors.Therefore, S stores domain-wall vectors representing a permutation if and only if it is an optimal solution that satisfies H nn a (X) = 2n.Table 3 presents the features of QUBO and Ising models derived from the expansion of mathematical expressions E nn a (X) and H nn a (S).The number of quadratic terms in these models is halved compared to conventional one-hot encoding.However, these models include linear terms with significantly large absolute values.The maximum absolute values for the QUBO and Ising models are 2n − 3 and n − 1, respectively.

QUBO/Ising Kernels with a Quadratic Number of Quadratic Terms Based on Our Dual-Matrix Domain-Wall Technique
Our new permutation encoding technique utilizes the inverse permutation.For a permutation π of n numbers, let π −1 denote the inverse, such that π −1 (π(i)) = i for all i (0 ≤ i ≤ n − 1).We refer to the permutations [π(0), π(1), . . ., π(n − 1)] and Recall that in conventional one-hot encoding, each row i (0 ≤ i ≤ n − 1) of an n × n matrix stores π(i) as a one-hot vector, as illustrated in Figure 1.It is easy to confirm that each column j (0 ≤ j ≤ n − 1) also stores a one-hot vector representing the inverse π −1 (j).Therefore, the row one-hot vectors and column one-hot vectors represent dual permutations.Figure 1 shows an example of a 4 × 4 matrix storing the row permutation [1, 3, 2, 0] and the column permutation [3, 0, 2, 1], which are dual permutations.Our dualmatrix domain-wall technique was inspired by this fact.It employs two matrices: matrix A = (a i,j ) (0 ≤ i ≤ n − 1 and 0 ≤ j ≤ n − 2) of size n × (n − 1) and matrix B = (b i,j ) (0 ≤ i ≤ n − 2 and 0 ≤ j ≤ n − 1) of size (n − 1) × n.These matrices are used to store dual permutations, where the rows of A represent a permutation and the columns of B represent its inverse permutation.For a clearer understanding, the reader should refer to Figure 4, which provides an example for n = 4 using a 4 × 3 matrix A and a 3 × 4 matrix B.  For domain-wall encoding, fixed guard bits/qubits are attached to the leftmost and rightmost columns of matrix A, as well as the top and bottom rows of matrix B, as illustrated in Figure 4. Specifically, we set a i,−1 = 1/ + 1 and a i,n−1 = 0/ − 1 for all i (0 ≤ i ≤ n − 1), and b −1,j = 1/ + 1 and b n−1,j = 0/ − 1 for all j (0 ≤ i ≤ n − 1) for QUBO/Ising models.Our goal is to design QUBO/Ising models that have optimal solutions satisfying the following three conditions: Row domain-wall condition: Each row of matrix A stores a domain-wall vector.
Column domain-wall condition: Each column of matrix B stores a domain-wall vector.
Dual-permutation condition: The row permutation of matrix A and the column permutation of matrix B are dual to each other.
An optimal solution of the designed QUBO/Ising models corresponds to any one of the possible n! permutations.
We begin by introducing the dual-matrix domain-wall technique for QUBO models.We define matrices ∆A = (∆a i,j ) and ∆B = (∆b i,j ) of size n × n (0 ≤ i, j ≤ n − 1) by computing the differences in each row of matrix A and each column of matrix B as follows: Figure 4 illustrates examples of ∆A and ∆B with n = 4.The QUBO kernel for generating a permutation using the dual-matrix domain-wall technique is defined by the following expression, denoted as E nn d (A, B): The first summation term takes the minimum value of 1 2 n if and only if the row domainwall condition is satisfied.Similarly, the second summation term takes the minimum value of 1  2 n if and only if the column domain-wall condition is satisfied.Lastly, the third summation term takes the minimum value of 0 if and only if ∆a i,j = ∆b i,j for all i and j.In essence, this condition signifies the satisfaction of the dual-permutation condition.It is important to note that these three summation terms can all reach their minimum values simultaneously.Consequently, the energy function E nn d (A, B) reaches its minimum value of n if and only if matrices A and B store dual permutations, satisfying all three conditions: the row domain-wall condition, the column domain-wall condition, and the dual-permutation condition.
For the Ising model, we can use the exact same mathematical expression as the QUBO model, resulting in the Hamiltonian function H nn d (A, B): Both the first and second summation terms take the minimum value of 2n if and only if both the row and column domain-wall conditions are satisfied.The third summation term takes the minimum value of 0 if and only if the dual-permutation condition is satisfied.Since these three summation terms can take the minimum values at the same time, the  (A, B) represents the QUBO model and takes the minimum value of n when all three conditions (row domain-wall, column domain-wall, and dual-permutation) are met.On the other hand, the Hamiltonian function H nn d (A, B) corresponds to the Ising model and achieves the minimum value of 4n when matrices A and B satisfy the dual-permutation condition along with the row and column domain-wall conditions.These QUBO/Ising models can be used as kernels for solving permutation-based combinational optimization problems as QUBO/Ising models.
The features of QUBO/Ising models obtained using our dual-matrix domain-wall encoding technique are presented in Table 3.It is evident from the table that these models consist of only 6n 2 − 12n + 4 quadratic terms.In contrast, models designed using the conventional one-hot encoding approach contain n 3 − n 2 quadratic terms, while models created using the all-different domain-wall encoding method have 1  2 n 3 − 3 2 n quadratic terms.Thus, our dual-matrix domain-wall technique significantly reduces the number of quadratic terms in the models.
Moreover, our encoding technique also leads to a reduction in the linear term coefficients in the Ising models.In Ising models obtained using conventional one-hot encoding and all-different domain-wall encoding, the linear terms have coefficients of 2n − 4 and n − 1, respectively.However, Ising models obtained through dual-matrix domain-wall encoding feature linear terms with a coefficient of 2.

QUBO/Ising Kernels for Partial Permutations
For two positive integers satisfying m ≤ n, a partial permutation refers to a sequence of m numbers selected from a set 0, 1, . . ., n − 1 of n numbers without repetition.A permutation is a special case of a partial permutation where m = n.In this section, we primarily focus on the case where m < n.However, we can also consider a permutation of n numbers by substituting m with n.A partial permutation can be represented by an injection π : {0, 1, . . ., m − 1} → {0, 1, . . ., n − 1}, where π(i) denotes the i-th element in the partial permutation [π(0), π(1), . . ., π(m − 1)].The total number of possible partial permutations is given by n!/(n − m)!.An inverse of a partial permutation π is a surjection p : {0, 1, . . ., n − 1} → {0, 1, . . ., m − 1} that satisfies p(π(i)) = i for all i (0 ≤ i ≤ m − 1).It should be noted that p(j) can take any value when π −1 (j) is undefined, meaning that there is no i satisfying π(i) = j.We refer to such π and p as dual permutations.
This section first demonstrates partial-permutation encoding using one-hot vectors.Subsequently, we illustrate the application of the dual-matrix domain-wall technique to generate partial permutations.It is important to note that the all-different domain-wall technique cannot be applied in this case.Unlike permutation generation, the column-wise sums of the domain-wall encoding for a partial permutation are not fixed.

Partial-Permutation Encoding by One-Hot Vectors
We can construct QUBO/Ising kernels to generate partial permutations using zeroone-hot vectors.For this purpose, we utilize an m × n matrix (m < n) of bits/qubits.Each row i (0 ≤ i ≤ m − 1) in the matrix represents π(i) as a one-hot vector, while each column j (0 ≤ j ≤ n − 1) represents π −1 (j) as a zero-one-hot vector.Here, we write π −1 (j) = φ if there is no i satisfying π(i) = j.Figure 5 illustrates an example of a 3 × 5 matrix that represents a partial permutation.In this example, the rows correspond to the partial permutation [π(0), π(1), π(2)] = [3, 0, 1] selected from the set {0, 1, 2, 3, 4}, and the columns represent the inverse [π −1 (0), (2) Ising model We can design the QUBO/Ising kernels E mn 1 (X)/H mn 1 (S) generating a partial permutation as follows: H mn 1 (S) = The first summation term in each expression takes the minimum value of 0 when every row is a one-hot vector.Similarly, the second summation term is evaluated as 0 when every column is a zero-one-hot vector.Consequently, the optimal solutions of E mn 1 (X)/H mn 1 (S) represent dual permutations.
We should note that previous works such as [26,27] have already presented QUBO models for partial permutation using zero-one-hot vectors.However, the mathematical expressions used in these works differ.In [26], a slack variable y i (0 ≤ i ≤ n − 1) was introduced for each column i, and the expression was utilized to ensure that each column represented a zero-one-hot vector.Similarly, in [27], the expres- was employed for the same purpose.In contrast, our proposed mathematical expression x i,j for zero-one-hot vectors is more intuitive and simpler as it does not yield any linear terms.

Partial-Permutation Encoding by Our Dual-Matrix Domain-Wall Encoding Technique
We utilize two matrices, namely A = (a i,j ) (0 ≤ i ≤ m − 1 and 0 ≤ j ≤ n − 2) of size m × (n − 1) with bit/qubit variables and B = (b i,j ) (0 ≤ i ≤ m − 1 and 0 ≤ j ≤ n − 1) of size (m − 1) × n.These matrices are employed for our dual-matrix domain-wall technique and represent partial permutations.For an illustration, refer to Figure 6, which shows example matrices A and B with dimensions m = 3 and n = 5, respectively.Similarly to dual-matrix domain-wall encoding for permutations, we include fixed guard bits/qubits in the leftmost and rightmost columns of A, as well as the top and bottom rows of B, as depicted in Figure 6.By applying Equation ( 16), we can derive two matrices ∆A = (∆a i,j ) and ∆B = (∆b i,j To satisfy the three conditions, we make modifications to the mathematical expressions E nn d (A, B) and H nn d (A, B) for the QUBO and Ising kernels that generate permutations.These modifications ensure that the row domain-wall condition, column domain-wall condition, and dual-permutation condition are met by the two matrices A and B with sizes m × (n − 1) and (m − 1) × n, respectively.The adjusted expressions are as follows: We proceed to demonstrate that E mn d (A, B) attains its minimum value of n − m if and only if ∆A represents a partial permutation.The first summation term takes the minimum value of 1  2 m when all rows of matrix A are domain-wall vectors, satisfying the row domainwall condition.Similarly, the second summation terms will reach the minimum values of 1 2 n if the column domain-wall condition is satisfied.Suppose that the row domain-wall condition and the column domain-wall condition are satisfied.Since the numbers of 1s in ∆A and ∆B are m and n, respectively, the third summation term takes the minimum value of 1  2 (n − m) if the dual-permutation condition is satisfied.Hence, E mn d (A, B) takes the value = n if the three conditions are satisfied.Unfortunately, this does not imply that E mn d (A, B) takes the minimum value of n if and only if the three conditions are satisfied.The third summation term can be smaller than 1  2 (n − m), say 0, if the row or column domain-wall condition is not satisfied.Hence, we need to prove that E mn d (A, B) is larger than n if at least one of the three conditions is not satisfied.More specifically, we prove the following lemma: Lemma 1. E mn d (A, B) is larger than n if A and B do not satisfy at least one of the row domain-wall, column domain-wall, and dual-permutation conditions.
Proof.Suppose that the row and column domain-wall conditions are satisfied and the numbers of 1s in ∆A and ∆B are m and n, respectively.We can observe that the third sum-mation term of E mn d (A, B) exceeds 1 2 (n − m) when ∆A and ∆B are not dual.Consequently, we can determine that E mn d (A, B) > 1 2 m + 1 2 n + 1 2 (n − m) = n.Next, let us consider the case in which the row and/or column domain-wall conditions are not fulfilled.If a row of A is not a domain-wall vector, the numbers of 1/−1s in ∆A are (1 + k)/k for some k ≥ 1.For example, if a row of A is 11001011, the corresponding row of ∆A would be 0010(−1)1(−1)100, containing three 1s and two −1s.Based on this observation, we can assign non-negative integers k A and k B such that the numbers of 1/−1s in ∆A and ∆B are (m + k A )/k A and (n + k B )/k B , respectively.It follows that ∆A and ∆B have m + 2k A and n + 2k B non-zero elements, respectively.Since we assume that the row and/or column domain-wall conditions are not satisfied, at least one of k A and k B must be greater than or equal to 1. Clearly, the first and second summation terms take values of 1 2 m + k A and 1  2 n + k B , respectively.The third summation term must be at least We can now analyze the value of E mn d (A, B) in the following cases: Case 1: The lower bound of E mn d (A, B) can be evaluated by the sum of the three summation terms as follows: Case 2: Similarly, E mn d (A, B) can be evaluated by the sum of the first and second summation terms as follows: Thus, this completes the proof of the lemma.
Since we have shown that E mn d (A, B) = n if the three conditions are satisfied, this lemma implies the following theorem: Theorem 1. E mn d (A, B) takes the minimum value of n if and only if A and B satisfy the row domain-wall, column domain-wall, and dual-permutation conditions are satisfied.
Next, we will evaluate the value of H mn d (A, B) when the three conditions are satisfied.The first summation term and second summation term take minimum values of 2m and 2n, respectively, if and only if the row and column domain-wall conditions are satisfied.Additionally, if the dual-permutation condition is satisfied, the third summation term takes a value of 2(n − m).Hence, H mn d (A, B) is equal to 2m + 2n + 2(n − m) = 4n if the row domain-wall, column domain-wall, and dual-permutation conditions are all satisfied.Similarly, we can prove that H mn d (A, B) is larger than 4n if at least one of the three conditions is not satisfied.Thus, we can state the following theorem: This theorem establishes the condition for H mn d (A, B) to attain its minimum value and confirms that it occurs when the row domain-wall, column domain-wall, and dualpermutation conditions are satisfied.
The readers should refer to Table 4, which provides the features of kernels E mn d (A, B) and H mn d (A, B).It can be observed that the number of quadratic terms in these models is given by 6mn − 4m − 4n.They utilize only a quadratic number of quadratic terms, whereas kernels E mn 1 (A, B) and H mn 1 (A, B) require a cubic number of quadratic terms.Additionally, the maximum absolute value of the linear term coefficients in H mn d (A, B) is only 2, while H mn 1 (A, B) requires large coefficients of m + n − 3.

Dual-Matrix Domain-Wall Technique Extended with a One-Hot Matrix
When the QUBO kernels designed by the dual-matrix domain-wall technique are involved in QUBO models for solving permutation-based combinatorial optimization problems, ∆A = (∆a i,j ) is used to compute the objective function of the optimization problems.For example, when checking if both π(i) = j and π(i ′ ) = j ′ are satisfied, a quadratic term ∆a i,j ∆a i ′ ,j ′ is employed.However, since ∆a i,j s are not variables, they are expanded using the variable a i,j s as shown below: Since a single quadratic term ∆a i,j ∆a i ′ ,j ′ is expanded to the four quadratic terms involving the variable a i,j s, the resulting QUBO model obtained through reduction contains numerous quadratic terms.To mitigate this issue, we introduce a matrix X = (x i,j ) with n × n-bit variables, which satisfies X = ∆A, into our QUBO kernel.Similarly, a one-hot matrix S = (s i,j ) with n × n-qubit variables is added to the Ising kernel.The matrix ensures that s i,j = ∆a i,j − 1 for all i and j, because each ∆a i,j takes 0 or 2. Since these matrices X and S store a permutation as a one-hot encoding, we refer to them as one-hot matrices To extend our dual-matrix domain-wall technique for one-hot matrices, we introduce an additional condition along with the row domain-wall condition, column domain-wall condition, and dual-permutation condition.This new condition is defined as follows: One-hot permutation condition.Each row of matrix X/S represents the permutation corresponding to matrix A in the form of a one-hot vector.
We begin by designing a QUBO model that satisfies these four conditions.The mathematical expression for this model is as follows: Similarly to the previous model, the first and second summation terms achieve the minimum value of n 2 when the row and column domain-wall conditions are satisfied, respectively.The third and fourth summation terms reach the minimum value of 0 if both the dual-permutation condition and the one-hot permutation condition are met.Therefore, matrix X stores a permutation if and only if E nn e (A, B, X) reaches its minimum value of n.
The same technique can also be applied to the Ising model.In this case, we introduce matrix S = (s i,j ) of size n × n, where s i,j + 1 = ∆a i,j holds for all i and j.The Ising model H nn e (A, B, S) can be designed as follows: Similarly to the QUBO model, the four conditions are satisfied, and matrix S stores a one-hot encoding of a permutation if and only if H nn e (A, B, S) reaches its minimum value of 4n.Optimal solutions of H nn e (A, B, S) correspond to a one-hot encoding stored in matrix S. The features of QUBO/Ising models obtained by expanding E nn e (A, B, X) and H nn e (A, B, S) are presented in Table 3.
We further apply the same technique for a partial permutation.However, a straightforward modification of E nn e (A, B, X)/H nn e (A, B, S) to fit them to a partial permutation may not generate QUBO/Ising models correctly.The optimal solutions of such QUBO/Ising models may not satisfy the four conditions.We modify their fourth summation term and obtain the following QUBO model E mn e (A, B, X): Note that the third and fourth summation terms do not have a coefficient of 1 2 , because the terms in the expanded formula will have a non-integer coefficient if they have a coefficient of 1  2 .The first and second summation terms take minimum values of 1 2 m and 1 2 n, respectively, if and only if the row and column domain-wall conditions are satisfied.The third summation term takes the minimum value of 0 if X and ∆A are equal.The fourth summation term takes the minimum value of 0 if x i,j = 0 or ∆b i,j = 1 holds for all i and j.
Thus, E mn e (A, B, X) = 1 2 m + 1 2 n if the four conditions are satisfied.We present the following lemma to prove that the "if" in the previous statement is "if and only if": Lemma 2. E mn e (A, B, X) is larger than 1 2 (m + n) if matrices A, B, and X do not satisfy at least one of the row domain-wall, column domain-wall, dual-permutation, and one-hot permutation conditions.
Proof.For any A, B, and X, the first, second, third, and fourth summation terms in E mn e (A, B, X) are at least 1 2 m, 1 2 n, 0, and 0, respectively.If the row/column domain-wall conditions are not satisfied, the first and second summation terms are larger than 1  2 m and 1 2 n, respectively.Hence, if the row and/or the column domain-wall conditions are not satisfied, E mn e (A, B, X) > 1 2 (m + n).Suppose that the row and column domain-wall conditions are satisfied.Clearly, matrices A and B have m and n 1s, respectively, and all the other elements are 0. Thus, the first and second summation terms take minimum values of 1  2 m and 1 2 n, respectively.If the dual-permutation condition is not satisfied, that is, ∆A and ∆B are not dual, then there exist i and j such that ∆a i,j = 1 and ∆b i,j = 0.For such an i and j, either (x i,j − ∆a i,j ) 2 or x i,j • (1 − ∆b i,j ) must be 1.Thus, at least one of the third or fourth summation terms is larger than 0. If the one-hot permutation condition is not satisfied, then X and ∆A are not equal, and the third summation term is larger than 0. Thus, E mn e (A, B, X) > 1 2 m + 1 2 n holds.This completes the proof of Lemma 2.
Since we have proved that E mn e (A, B) = 1 2 m + 1 2 n if the four conditions are satisfied, we can state the following theorem: Theorem 3. E mn e (A, B, X) takes the minimum value of 1 2 m + 1 2 n if and only if A and B satisfy the row domain-wall, column domain-wall, dual-permutation, and one-hot permutation conditions.
We can extend the Ising model for a partial permutation to have a one-hot matrix S in the same way as follows: If the four conditions are satisfied, then the first, second, third, and fourth summation terms take the values 2m, 2n, 0, and 0, respectively.Thus, H mn e (A, B, S) = 2m + 2n holds when the four conditions are satisfied.Similarly, we can prove that H mn e (A, B) is larger than 2m + 2n if at least one of the four conditions is not satisfied.Thus, we can state the following theorem: The reader should refer to Table 4, which provides the features of E mn e (A, B, X) and H mn e (A, B, S).It can be observed that the number of quadratic terms in these kernels is given by 6mn − 4m − 4n, while E mn d (A, B) and H mn d (A, B) have 6mn − 6m − 6n + 4 quadratic terms.Therefore, the introduction of one-hot matrices results in a relatively small increment in the quadratic term counts.This reduction in the number of quadratic terms highlights the effectiveness of utilizing one-hot matrices X and S in the QUBO and Ising models for partial permutations.By employing these techniques, we can significantly enhance the efficiency of solving optimization problems associated with partial permutations

Applications of QUBO/Ising Kernels for Generating Permutations in Combinatorial Optimization Problems
In this section, we begin by introducing the Particle Placement Problem (PPP), which serves as a fundamental problem that many permutation-based combinatorial optimization problems can be reduced to with ease.We demonstrate that the PPP can be effectively reduced to QUBO/Ising models with kernels for generating permutations.Furthermore, we show that several permutation-based combinatorial optimization problems, such as the quadratic assignment problem (QAP), traveling salesman problem (TSP), and the subgraph isomorphism problem, can be equivalently reduced to the PPP.Consequently, these combinatorial optimization problems can be reduced to QUBO/Ising models through the PPP.

The Particle Placement Problem (PPP)
Suppose we have m particles with IDs 0, 1, . . ., m − 1 and n positions with IDs 0, 1, . . ., n − 1, where m ≤ n is satisfied.Let us consider the PPP, which involves placing m particles in n positions without conflicts.The particle placement is determined by injection π : {0, 1, . . ., m − 1} → {0, 1, . . ., n − 1}, where each particle i (0 ≤ i ≤ m − 1) is placed in position π(i).The problem aims to find an optimal placement π that minimizes the total energy, which is defined by two types of energies as follows: Potential P i,j : the energy that occurs when particle i is placed in position j (0 ≤ i ≤ m − 1 and 0 ≤ j ≤ n − 1).
Interaction I i,j,i ′ ,j ′ : the energy that occurs between two particles i and i ′ when they are placed in positions j and j ′ , respectively, where (i, Here, Q(m, n) is a set of all consistent quartets (i, j, i ′ , j ′ ) satisfying 0 ≤ i < i ′ ≤ m − 1, 0 ≤ j, j ′ ≤ m − 1, and j ̸ = j ′ .Both potential P i,j and interaction I i,j,i ′ ,j ′ can be any integer, including negative values.The total energy of π is the sum of all potentials and interactions, given by the equation The particle placement problem (PPP) aims to find an optimal placement π with the minimum PPP(π) among all possible particle placements.
We will demonstrate that the total energy of π can be computed using the corresponding binary one-hot matrix, and we can reduce the PPP to a QUBO problem.In this sub-section, we focus on the PPP with m particles and n locations, where m < n and π represents a partial permutation.However, it is worth noting that the same technique can be applied to the PPP with n particles and n locations, as it is relatively straightforward.
Suppose that an m × n-bit matrix X stores a partial permutation π as one-hot vectors.In other words, x i,j = 1 holds if and only if π(i) = j.Let PPP(X) be a quadratic formula defined as follows: Clearly, x i,j = 1 if and only if π(i) = j, and x i,j x i ′ ,j ′ = 1 if and only if π(i) = j and π(i ′ ) = j ′ .Thus, we have PPP(X) = PPP(π).Since PPP(X) is a quadratic formula of X, we can design a QUBO model using a QUBO kernel E mn 1 (X) corresponding to the PPP as follows: where λ is a parameter large enough to prioritize E mn 1 (X) and guarantee that X is a one-hot vector.By finding an n × n-bit matrix X that minimizes E PPP 1 (X), we can obtain the optimal solution π of the PPP corresponding to such an X.
We can use the other QUBO kernels E mn d (A, B) and E mn e (A, B, X) instead of E mn 1 (X): Similarly, by finding their optimal solutions, we can obtain the optimal solution π of the PPP.
For an m × n-qubit matrix S = (s i,j ), let S ′ = (s ′ i,j ) denote the matrix of the same size satisfying s ′ i,j = s i,j + 1 for all i and j.Clearly, s ′ i,j = 0/2 if s i,j = −1/+1, respectively.Thus, for PPP(S ′ ) computed by the mathematical expression in Equation (32), PPP(S ′ ) = 4PPP(π) is satisfied if S stores permutation π as one-hot vectors.Using this fact, we can design Ising models for solving the PPP as follows: H PPP e (A, B, S) = λH mn e (A, B, S) + PPP(S ′ ) It should be clear that the optimal solutions of these Ising models give the optimal solutions of the PPP if λ is large enough to prioritize the Ising kernels for generating partial permutations.
The PPP can have a total of mn potentials and 1 2 m 2 n 2 − 1 2 m 2 n − 1 2 mn 2 + 1 2 mn interactions.When the PPP is converted to the equivalent QUBO/Ising models, non-zero interactions may lead to quadratic terms.In essence, the total number of quadratic terms in the QUBO/Ising models for solving the PPP is the sum of the number of terms produced by a kernel generating permutations and the number of quadratic terms for non-zero interactions in the PPP.As demonstrated in this paper, our dual-matrix domain-wall technique significantly reduces the count of quadratic terms.However, if the number of quadratic terms derived from non-zero interactions in the PPP is much larger than that for the kernel generating permutations, the advantage of our dual-matrix domain-wall technique becomes limited.It is worth noting that both the conventional one-hot technique and the all-different domain-wall technique require a cubic number of quadratic terms, whereas our dual-matrix domain-wall technique only uses a quadratic number of quadratic terms.If the PPP has full non-zero interactions, it would result in a number of quadratic terms that is proportional to the fourth power.In such a scenario, the resulting QUBO/Ising models would have a fourth power number of quadratic terms, and the reduction effect achieved by the dual-matrix domain-wall technique for the kernel terms would be relatively small.To evaluate the effect of the dual-matrix domain-wall technique on specific combinatorial optimization problems, we will present several examples and assess the number of non-zero interactions in the PPP.

Combinatorial Optimization Problems that Can Be Reduced to the PPP
This subsection demonstrates that several combinatorial optimization problems, namely the quadratic assignment problem (QAP), the traveling salesman problem (TSP), and the sub-graph isomorphism problem, which are known to be NP-hard [33], can be effectively reduced to the PPP.Consequently, the instances of these problems can be transformed into equivalent QUBO/Ising models.Additionally, we illustrate that the maximum weight matching problem can also be reduced to the PPP.While this problem is not NP-hard and can be solved in polynomial time [34], there is currently no known efficient parallel algorithm for general graphs [35].Consequently, if a highly powerful QUBO/Ising solver becomes available in the future, it could be advantageous to solve this problem by utilizing the reduction to QUBO/Ising models.

The Quadratic Assignment Problem (QAP)
The quadratic assignment problem (QAP) [36] aims to find the optimal arrangement π of n facilities to n positions.This problem involves two types of instance parameters: flows and distances.Let f i,i ′ and d j,j ′ denote the flow from facility i to i ′ (0 ≤ i, i ′ ≤ n − 1) and the distance from j to j ′ (0 ≤ j, j ′ ≤ n − 1), respectively.The objective of QAP is to find a permutation π that minimizes the total logistics cost, defined by the following formula: To reduce this problem to the PPP, we set ,j for all i, i ′ , j, and j ′ ((i, j, i ′ , j ′ ) ∈ Q(n, n)) and set P i,j = f i,i • d j,j for all i and j (0 ≤ i, j ≤ n − 1).Since QAP(π) = PPP(π) holds for all permutation π, any QAP instance can be reduced to the PPP.The reduced PPP has up to 1  2 n 4 − n 3 + 1 2 n 2 non-zero interactions and n 2 non-zero potentials.If not all flows are non-zero, the number of non-zero interactions may be reduced.Let k be the number of pairs of facilities i and i ′ (0 ≤ i < i ′ ≤ n − 1) such that at least one of f i,i ′ or f i ′ ,i is non-zero.Clearly, k is at most 1 2 n 2 − 1 2 n.In this case, the reduced PPP has only 1  2 n 2 k − 1 2 nk non-zero interactions.

The Traveling Salesman Problem (TSP)
The traveling salesman problem (TSP), which aims to find the shortest route to visit all n cities with IDs 0, 1, . .., n − 1 can be reduced to the PPP as follows.Let d j,j ′ (0 ≤ j, j ′ ≤ n − 1) denote the distance between two cities j and j ′ .We use a permutation π for representing a route, where the i-th visited city (0 ≤ i ≤ n − 1) corresponds to the city with ID π(i).The objective of the TSP is to find a permutation π that minimizes the total route length computed using the following formula: To reduce the TSP to the PPP, we define I i,j,i ′ ,j ′ = d j,j ′ for all i, i ′ , j, and j ′ if i ′ = (i + 1) mod n.For all other I i,j,i ′ ,j ′ and P i,j , we set them to 0. This reduction allows us to observe that TSP(π) = PPP(π) holds for all permutations π.Therefore, it is possible to reduce any TSP instance to the PPP.The reduced PPP has n 3 − n 2 non-zero interactions.
We can consider the traveling salesman problem (TSP) for a weighted undirected graph in which each edge is assigned a weight value.The goal of the TSP is to find a Hamiltonian cycle in the graph with the minimum total weight that visits all nodes.Figure 7 shows an example of a 40-node weighted graph.We assume that the weight of each edge is its length in the figure.The figure also depicts the optimal solution of the TSP.Note that if an input graph has no Hamiltonian cycle, the TSP has no feasible solution.Suppose that we have a weighted graph with nodes 0, 1, . .., n − 1, and let d j,j ′ (0 ≤ j, j ′ ≤ n − 1) be the weight of an edge (j, j ′ ).If the graph does not have an edge (j, j ′ ), we assume d j,j ′ = +∞ to prevent selecting (j, j ′ ) as an edge of a cycle.The total route of a node permutation π is computed using Equation (40), based on the defined d j,j ′ values.For the real implementation, we set d j,j ′ = BIG instead of +∞, where BIG is sufficiently larger than the maximum edge weight.An edge (j, j ′ ) with d j,j ′ = BIG is selected only if the TSP has no feasible solution.
We will introduce the technique for reducing the TSP into the PPP, so that the reduced PPP becomes sparse.For this purpose, we apply a fixed bias of −BIG, so that I i,j,i ′ ,j ′ = d j,j ′ − BIG for all i, i ′ , j, and j ′ if i ′ = (i + 1) mod n.For all other I i,j,i ′ ,j ′ and P i,j , we set them to 0. This reduction allows us to observe that TSP(π) = PPP(π) + n • BIG holds for all permutations π.Additionally, I i,j,i ′ ,j ′ is non-zero only if the graph has an edge (j, j ′ ) and i ′ = (i + 1) mod n holds.Thus, the PPP becomes so sparse that it has only 2en non-zero interactions, where e is the number of edges in the graph.For instance, if no two edges intersect each other, as illustrated in Figure 7, the graph has at most 3n − 6 edges.Hence, the reduced PPP has at most 2(3n − 6)n = 6n 2 − 12 interactions, which is much smaller than the original TSP with n 3 − n 2 interactions.

The Sub-graph Isomorphism Problem
Let G = (V, E) and G ′ = (V ′ , E ′ ) be guest and host graphs such that V = {0, 1, . . ., m − 1} and V ′ = {0, 1, . . ., n − 1}.The sub-graph isomorphism problem aims to find a permutation π such that (π(i), π(i ′ )) ∈ E ′ holds for all (i, i ′ ) ∈ E. If such a permutation π exists, it implies that G is a sub-graph of G ′ .To reduce the sub-graph isomorphism problem to the PPP, we assign −1 to I i,j,i ′ ,j ′ and I i,j ′ ,i ′ ,j if both (i, i ′ ) ∈ E and (j, j ′ ) ∈ E ′ hold, and assign 0 to all other I i,j,i ′ ,j ′ as well as all P i,j .It is evident that I i,π(i),i ′ ,π(i ′ ) = −1 holds for all edges (i, i ′ ) in G if π is a permutation that proves G is a sub-graph of G ′ .Therefore, PPP(π) = −|E| if such a permutation exists.Consequently, any instance of the sub-graph isomorphism problem can be reduced to the PPP.The reduced PPP has 2 Let G = (V, E, w) be a weighted graph, where V = {0, 1, . . ., n − 1} and w i,j represents the weight of an edge (i, j) (i < j) in E. We assume that w i,j = 0 for all other pairs (i, j) such that i ≥ j or (i, j) is not in E. A subset of E is called a matching if none of its edges share a common node.The maximum weight matching problem aims to find a matching M (⊆ E) with the largest total weight.We can represent a matching M using a permutation π as follows.An edge (i, j) (∈ E) is in the matching M if both π(i) = j and π(j) = i hold, and node i is not connected to any edge in M if π(i) = i.To transform the maximum weight matching problem into the PPP, we assign −w i,j to I i,j,j,i for all (i, j) ∈ E, while assigning 0 to all other I i,j,i ′ ,j ′ as well as all P i,j .Let π be a permutation that represents a matching M. Since w i,π(i) = 0 whenever i ≥ π(i), we can compute the total weight of M using the obtained PPP as follows: Therefore, any instance of the maximum weight matching problem can be reduced to the PPP with |E| interactions.If G is a connected graph, n − 1 ≤ |E| ≤ 1 2 n 2 − 1 2 n holds, and so the number of interactions is from n − 1 to 1  2 n 2 − 1 2 n.

The Bipartite Maximum Weight Matching Problem
If an instance of the maximum weight matching problem is restricted to a bipartite graph, a more efficient reduction to the PPP can be used.Consider a weighted bipartite graph, denoted as G = (V, V ′ , E, w), where V = {0, 1, . . ., m − 1}, V ′ = {0, 1, . . ., n − 1}, E ⊆ V × V ′ , and w i,j represents the weight of an edge (i, j) ∈ E. We can use a partial permutation π : V → V ′ to represent a matching M of G such that (i, π(i)) ∈ M. To transform the maximum weight matching problem for the weighted bipartite graph into the PPP, we assign −w i,j to P i,j for all (i, j) ∈ E, while assigning 0 to all other P i,j as well as all I i,j,i ′ ,j ′ .Thus, we can compute the total weight of M using the obtained PPP as follows: Therefore, any instance of the bipartite maximum weight matching problem can be reduced to the PPP.The resulting PPP has |E| non-zero potentials but has no non-zero interactions.

The Quadratic Term Counts of the Ising Model
This section presents the quadratic term counts for various permutation-based combinatorial optimization problems.We examined the total number of quadratic terms in Ising models and permutation kernels to observe the impact of smaller permutation kernels on the size of Ising models.To facilitate this analysis, we introduced the concept of PPP density, which represents the ratio of non-zero interactions to the maximum number of interactions |Q(m, n)| = 1 2 m 2 n 2 − 1 2 m 2 n − 1 2 mn 2 + 1 2 mn.We used PyQUBO [37], a Python tool designed to generate QUBO/Ising models from mathematical expressions, and SymPy [32], a Python library for symbolic mathematics, to validate the accuracy of the QUBO/Ising models.However, this approach, which operates on provided mathematical expressions, has a significant drawback in terms of its substantial time and space requirements.This limitation could make it impractical to solve particularly large problems.For optimization problems based on specific permutations, an alternative and more efficient approach is to create a custom computer program for generating QUBO/Ising models.This approach can leverage a low-level programming language like C++.Given the fixed and straightforward nature of the underlying mathematical expressions for such problems, it becomes feasible to design a program that efficiently generates QUBO/Ising models without the need to expand the mathematical expressions.This process can be achieved with linear time complexity and requires a memory proportional to the size of the models.
Figure 8 illustrates the quadratic term counts for the following problem instances: 1.
A random instance of the QAP with n = 50: The reduced PPP consists of all 3,001,250 interactions, yielding a density of 1.00.

2.
A random instance of the QAP with n = 50 and 10% non-zero flows: The reduced PPP contains 301,350 interactions out of 3,001,250 maximum interactions, yielding a density of 0.10.

3.
A random instance of the TSP with n = 100: The reduced PPP includes 990,000 interactions out of 49,005,000 maximum interactions, yielding a density of 0.020.

4.
The TSP for a randomly generated 40-node weighted graph as shown in Figure 7: The graph has 104 edges, and the reduced PPP includes 8320 interactions out of 1,216,800 maximum interactions, yielding a density of 0.0068.

5.
The TSP for a randomly generated 300-node weighted graph: The graph is generated to avoid edge intersections, similar to the 40-node weighted graph.It has 277 edges, and the reduced PPP includes 524,400 interactions out of a maximum of 4,023,045,000 interactions, yielding a density of 0.00013.

6.
A random instance of the sub-graph isomorphism problem for a three-regular guest graph with m = 200 nodes and a six-regular host graph with n = 400 nodes: The reduced PPP consists of 720,000 interactions out of 3,176,040,000 maximum interactions, yielding a density of 0.00023.7.
A random instance of the maximum weight matching problem for a complete graph with n = 300 nodes: The reduced PPP contains 44,850 interactions out of 4,023,045,000 maximum interactions, yielding a density of 0.000011.8.
A random instance of the bipartite maximum weight matching problem for a complete bipartite graph with m = n = 300 nodes: The reduced PPP has no interactions, yielding a density of zero.
The figure illustrates the total number of quadratic terms in the resulting Ising model, referred to as the "model total", for each problem.Additionally, it displays the number of quadratic terms contributed by the permutation kernels, referred to as the "permutation kernel", which is also included in the model total.The quadratic term counts are evaluated for four permutation kernels: the conventional one-hot technique (one-hot), the all-different domain-wall technique (all-different), the dual-matrix domain-wall technique (dual-matrix), and the extended technique with a one-hot matrix (extended).problems by the conventional one-hot technique (one-hot), the all-different domain-wall technique (all-different), the dual-matrix domain-wall technique (dual-matrix), and that extended with a one-hot matrix (extended).We evaluated the total quadratic term counts of the resulting Ising models (model total), as well as those specifically for kernels generating permutations (permutation kernel).
From the figure, it is evident that our dual-matrix domain-wall technique is more effective in reducing the quadratic terms of resulting Ising models when the original PPPs have a lower density.When the density of the PPP is close to 1, the quadratic terms for PPP interactions outweigh those for permutation kernels.On the other hand, when the PPP has a lower density, the total number of quadratic terms in the resulting Ising model can be significantly reduced.For instance, if the conventional one-hot encoding is used for the TSP with a randomly generated 300-node graph, the total number of quadratic terms in the resulting Ising model is 27,434,400.However, when our dual-matrix domain-wall technique extended by a one-hot matrix is applied, the total number of quadratic terms reduces to only 1,062,000.This represents a remarkable reduction factor of 25.8.
We should take note of the impact of extending by a one-hot matrix in our dual-matrix domain-wall technique.Recall that the original technique utilizes variables in matrix A for PPP interactions, while the extended technique introduces an one-hot matrix X/S.As shown in Equation ( 26), the dual-matrix domain-wall technique generates four quadratic terms of A for each PPP interaction.In contrast, the extended dual-matrix domain-wall technique produces only one quadratic term.When the PPP density is high and involves numerous interactions, the quadratic terms generated by A may overlap, resulting in the Ising model encompassing nearly all pairs of variables in A. Similarly, the Ising models produced by the extended dual-matrix domain-wall technique also include almost all pairs of variables in X/S.Consequently, the number of quadratic terms is nearly identical for both techniques, and the extended dual-matrix domain-wall technique does not effectively reduce them.Likewise, when the PPP density is too low and contains few interactions, the number of quadratic terms generated by variables A/X/S is small.Since the extended dualmatrix domain-wall technique yields a slightly larger kernel size, it does not perform well in reducing the overall number of quadratic terms.Therefore, the extended dual-matrix domain-wall technique can effectively reduce the total number of quadratic terms only when the PPP density is moderate-neither too small nor too large.An example illustrating this is the sub-graph isomorphism problem depicted in Figure 8, where the extended dual-matrix domain-wall technique reduces the total number of quadratic terms to less than half.The quadratic term counts of the extended technique are slightly larger than those of the non-extended technique for the other problem.However, since the difference is quite small, we should select the extended dual-matrix domain-wall technique if we must use a fixed technique for some reason.

Embedding of Permutation Kernels in D-Wave Quantum Annealer
This section presents the implementation results of Ising kernels on a quantum annealer called the D-Wave Advantage 4.1.The purpose was to evaluate the suitability of the dual-matrix domain-wall technique for currently available quantum annealers and to verify the impact of quadratic term counts in Ising models on the cell counts.We employed Ising kernels designed using mathematical expressions such as H nn 1 (S) (one-hot), H nn a (S) (all-different), H nn d (S) (domain-wall), and H nn e (S) (extended) for generating permutations.Additionally, we used H mn 1 (S) (one-hot), H mn d (S) (domain-wall), and H mn e (S) (extended) for partial permutations.These Ising kernels were embedded in the D-Wave Advantage 4.1, which was equipped with a 5760-node Pegasus graph designed for solving Ising models.However, due to faulty nodes and connections, only 5627 cells with 40,279 connections were available.Therefore, we retrieved the real graph topology and employed "minorminer", a heuristic tool for finding minor embeddings [38], to embed the Ising kernels.
Figure 9 illustrates the number of cells used to embed Ising kernels in the D-Wave Advantage 4.1.We embedded Ising kernels to generate permutations obtained by H nn 1 (S) (one-hot), H nn a (S) (all-different), H nn e (S) (domain-wall), and H nn d (S) (extended) for n = 5 until embedding failed.For partial permutations, we fixed the parameter m = 12 and evaluated embedding for H mn 1 (S) (one-hot), H mn d (S) (domain-wall), and H mn e (S) (extended) for n = 12 until embedding failed.Ising kernels for generating permutations using the one-hot, all-different, domain-wall, and extended techniques were successfully embedded in the D-Wave Advantage 4.1 topology until n = 16, 28, 35, and 29, respectively.Regarding Ising kernels for partial-permutation generation by the one-hot, domain-wall, and extended techniques, successful embedding was achieved until n = 22, 112, and 75, respectively.
Our dual-matrix domain-wall technique significantly reduced the cell count required for permutation generation.However, the permutation kernels produced by the extended dual-matrix domain-wall technique utilize slightly more cells than those produced by the non-extended version.Additionally, the total number of cells used by the extended technique can be smaller, because it uses fewer quadratic terms for PPP interactions.
If the entire Ising models for solving combinatorial optimization problems are embedded, the upper bound for the value of n must be smaller.For such small values of n, conventional digital computers are capable of finding optimal solutions for the problems.Therefore, it does not make sense to utilize the D-Wave Advantage 4.1 for solving such problems.However, in the future, if a larger number of cells becomes available, the D-Wave Advantage 4.1 could potentially load Ising models reduced from larger combinatorial optimization problems.In this scenario, it could uncover optimal solutions that are beyond the reach of conventional digital computers.In such cases, kernels for generating permutations with larger values of m and n would be employed, and our dual-matrix domain-wall technique would be a powerful approach for designing them.For a quantum annealer with a fixed size, our dual-matrix domain-wall technique enables us to program larger-sized combinatorial optimization problems than the conventional one-hot encoding techniques, which require a cubic number of quadratic terms.

Conclusions
We introduced a novel permutation encoding technique called the dual-matrix domain wall, which can be utilized to design QUBO/Ising kernels for generating permutations and partial permutations.The conventional one-hot encoding technique for QUBO/Ising kernels requires a cubic number of quadratic terms, whereas our dual-matrix domain-wall technique enables a reduction to a quadratic number of quadratic terms.In this study, we provided a detailed mathematical analysis of QUBO/Ising kernels for permutation generation.
Additionally, we introduced a generic problem known as the particle placement problem (PPP), which can be easily converted to QUBO/Ising models by utilizing the QUBO/Ising kernels for generating permutations.Furthermore, we identified several permutation-based combinatorial optimization problems that can be reduced to the PPP.As a result, these problems can be efficiently converted to QUBO/Ising models using our dual-matrix domain-wall technique.
To evaluate the effectiveness of our approach, we analyzed the number of quadratic terms used in the QUBO/Ising models derived from these optimization problems.Finally, we provided implementation results on the D-Wave Advantage 4.1, a quantum annealer developed by D-Wave Systems.
Regrettably, the existing quantum annealers currently accessible are of insufficient size to showcase the practical efficacy of our novel technique.From a practical perspective, our methodology is specifically tailored for optimal performance exclusively on quantum annealers of a substantial scale.The dual-matrix domain-wall technique is not presently feasible with the D-Wave Advantage 4.1.Nonetheless, it possesses the potential to evolve into a valuable technique for forthcoming large-scale quantum annealers.

1 Figure 3 .
Figure 3. Example of minor embedding for the Ising model of Figure 2 to the quantum annealer.

Figure 5 .
Figure 5. 3 × 5 matrix representing a partial permutation selecting 3 from 5 using one-hot and zero-one-hot vectors.

Theorem 2 .
H mn d (A, B) takes the minimum value of 4n if and only if A and B satisfy the row domain-wall, column domain-wall, and dual-permutation conditions.

Theorem 4 .
H mn e (A, B, S) takes the minimum value of 2m + 2n if and only if A, B, and S satisfy the row domain-wall, column domain-wall, dual-permutation, and one-hot permutation conditions.

Figure 7 .
Figure 7.A 40-node weighted graph with the optimal TSP solution.

Figure 8 .
Figure8.The quadratic term counts of Ising models reduced from combinatorial optimization problems by the conventional one-hot technique (one-hot), the all-different domain-wall technique (all-different), the dual-matrix domain-wall technique (dual-matrix), and that extended with a one-hot matrix (extended).We evaluated the total quadratic term counts of the resulting Ising models (model total), as well as those specifically for kernels generating permutations (permutation kernel).

Figure 9 .
Figure 9.The numbers of cells of the D-Wave Advantage 4.1 used for embedding permutation kernels.We evaluated them for Ising kernels that generate permutations of n numbers and partial permutations of m = 12 numbers selected from n numbers.
Hamiltonian function H nn d (A, B) takes the minimum value of 4n if and only if matrices A and B store dual permutations.It is worth noting that although E nn d (A, B) and H nn d (A, B) have the same mathematical expressions, their expansions differ.The energy function E nn d

Table 4 .
QUBO/Ising kernels for generating a partial permutation selecting m numbers from n numbers.