JP7795103B2 - Program, data processing device and data processing method - Google Patents

Description

The present invention relates to a program, a data processing device, and a data processing method.

Assignment problems are a class of NP-hard combinatorial optimization problems with a variety of real-world applications, such as vehicle routing and FPGA (Field-Programmable Gate Array) block placement. An example of an assignment problem is the quadratic assignment problem (QAP) (see, for example, Non-Patent Document 1).

The quadratic assignment problem is a problem of finding an assignment that minimizes the sum of the products of the flow volume between elements (costs such as transporting materials between facilities) and the distances between the destinations to which each element is assigned when allocating n elements (facilities, etc.) to n destinations. In other words, the quadratic assignment problem is a problem of finding an assignment that minimizes the value of the evaluation function expressed in equation (1) below. The evaluation function represents the cost according to the allocation state and is also called a cost function.

In equation (1), f _i,j indicates the flow amount between elements with identification numbers = i and j, and d _φ(i),φ(j) indicates the distance between the assignment destinations to which the elements with identification numbers = i and j are assigned.
Meanwhile, there is a Boltzmann machine (also called an Ising machine) that uses an Ising evaluation function and is a type of recurrent neural network that can solve large-scale discrete optimization problems, which von Neumann computers are not good at.

A Boltzmann machine converts a combinatorial optimization problem into an Ising model that represents the spin behavior of a magnetic material. The Boltzmann machine then searches for the state of the Ising model that minimizes the value of an Ising-type evaluation function (see, for example, Non-Patent Document 3) using Markov chain Monte Carlo methods such as simulated annealing and parallel tempering (see, for example, Non-Patent Document 2). The value of the evaluation function corresponds to energy. Note that by changing the sign of the evaluation function, the Boltzmann machine can also search for the state where the value of the evaluation function is maximized. The state of the Ising model can be expressed by a combination of the values of multiple state variables (also called neuron values). The value of each state variable can be either 0 or 1.

The Ising-type evaluation function is defined, for example, by the following equation (2):

The first item on the right-hand side is the sum of the products of the values (0 or 1) of two state variables selected from the N state variables, without omission or duplication, and the weight value (representing the strength of interaction between the two state variables), for all combinations of all state variables of the Ising model. s _i is the state variable with identification number i, s _j is the state variable with identification number j, and wi _,j is a weight value indicating the magnitude of interaction between the state variables with identification numbers i and j. The two items on the right-hand side are the sum of the products of the bias coefficient and the state variable for each identification number. _{b i} indicates the bias coefficient for identification number = i.

The amount of change in energy (ΔE _i ) associated with a change in the value of s _i is expressed by the following equation (3).

In equation (3), when s _i changes from 1 to 0, Δs _i becomes -1, and when s _i changes from 0 to 1, Δs _i becomes 1. Note that h _i is called a local field, and ΔE _i is obtained by multiplying h _i by the sign (+1 or -1) according to Δs _i .

Then, for example, if ΔE _i is smaller than a noise value (also called thermal noise) obtained based on the random number and the temperature parameter value, the value of s _i is updated to cause a state transition, and the local field is also updated, and this process is repeated.

Techniques have been proposed for solving quadratic assignment problems using such Boltzmann machines (see, for example, Patent Document 1 and Non-Patent Document 4).
The Ising-type evaluation function for the quadratic assignment problem can be expressed by the following equation (4).

In equation (4), x is a vector of state variables and represents the allocation state of n elements to n allocation destinations. ^xT can be expressed as ( _x1,1 , ..., x1 _,n , _x2,1 , ..., _x2,n , ..., xn _,1 , ..., xn _,n ). _xi,j = 1 indicates that the element with identification number = i is allocated to the allocation destination with identification number = j, and xi _,j = 0 indicates that the element with identification number = i is not allocated to the allocation destination with identification number = j.

W is a matrix of weight values, and can be expressed by the following equation (5) using the flow amount (f _i,j ) and the matrix D of the distances between n destinations.

US Patent Application Publication No. 2021/0326679

Eugene L. Lawler, “The quadratic assignment problem”, Management Science, Vol.9, No.4 pp.586-599, July 1963 Robert H. Swendsen and Jian-Sheng Wang, “Replica monte carlo simulation of spin-glasses”, Physical Review Letters, Vol.57, No.21, pp.2607-2609, November 1986 K. Dabiri, et al., “Replica Exchange MCMC Hardware With Automatic Temperature Selection and Parallel Trial”, IEEE Transactions on Parallel and Distributed Systems, Vol.31, No.7, pp.1681-1692, July 2020 M. Bagherbeik et al., “A permutational boltzmann machine with parallel tempering for solving combinatorial optimization problems”, In International Conference on Parallel Problem Solving from Nature, pp.317-331, Springer, 2020 Rainer E Burkard, Stefan E Karisch, and Franz Rendl, “Qaplib-a quadratic assignment problem library”, Journal of Global optimization, 10(4), pp.391-403, 1997 Gintaras Palubeckis, “An algorithm for construction of test cases for the quadratic assignment problem”, Informatica, Lith. Acad. Sci., Vol.11, No.3, pp.281-296, 2000 Zvi Drezner, Peter M Hahn, and Eric D Taillard, “Recent advances for the quadratic assignment problem with special emphasis on instances that are difficult for meta-heuristic methods”, Annals of Operations research, 139(1), pp.65-94, 2005 Allyson Silva, Leandro C. Coelho, and Maryam Darvish, “Quadratic assignment problem variants: A survey and an effective parallel memetic iterated tabu search”, European Journal of Operational Research, 2020 Danny Munera, Daniel Diaz, and Salvador Abreu, “Hybridization as cooperative parallelism for the quadratic assignment problem”, In Hybrid Metaheuristics, pp. 47-61, Springer International Publishing Switzerland, 2016 Kresimir Mihic, Kevin Ryan, and Alan Wood, “Randomized decomposition solver with the quadratic assignment problem as a case study”, INFORMS Journal on Computing, Vol.30, No.2, pp.295-308, 2018

Methods of calculating quadratic assignment problems using evaluation functions like those above require repeated calculations of energy changes and local field updates, as well as repeated memory accesses, which takes a long time to calculate.

In one aspect, the present invention aims to provide a program, data processing device, and data processing method that can quickly calculate assignment problems.

In one embodiment, a program is provided that causes a computer to execute a process of searching for a solution to an assignment problem through local search using an evaluation function that represents costs according to the assignment state. The program calculates, using vector arithmetic operations, a first change in the evaluation function when a first assignment change occurs in which the assignment destinations of a first element and a second element among the multiple elements are swapped, based on a flow matrix that represents the flow volume between multiple elements assigned to multiple assignment destinations and a distance matrix that represents the distances between the multiple assignment destinations, both stored in memory; determines whether to allow the first assignment change based on the first change; and, if it is determined that the first assignment change is allowable, updates the assignment state and updates the distance matrix so that two columns or two rows corresponding to the first element and the second element are swapped.

Also, in one embodiment, a data processing apparatus is provided.
Also, in one embodiment, a data processing method is provided.

In one aspect, the present invention allows for fast computation of assignment problems.

FIG. 10 is a diagram illustrating an example of calculating QAP. 10A and 10B are diagrams illustrating an example of rearrangement of a distance matrix when calculating a QAP and an example of a data processing device. FIG. 10 is a diagram illustrating an example of update calculation of a cache matrix. FIG. 1 illustrates an example of a solver system that performs parallel tempering. FIG. 10 is a diagram illustrating an example of an algorithm for searching for a solution to a QAP by local search using parallel tempering. 10 is a flowchart showing the overall processing flow of local search by parallel tempering. 10 is a flowchart illustrating an example of a process for initializing a replica in the case of a QAP. 10 is a flowchart illustrating an example of a replica search process in the case of a QAP. FIG. 10 is a diagram illustrating an example of an algorithm for searching for a solution to QSAP by local search using parallel tempering. 10 is a flowchart illustrating an example of the flow of a replica initialization process in the case of QSAP. 10 is a flowchart illustrating an example of a replica search process in the case of QSAP. FIG. 10 is a diagram showing the evaluation results of the degree of speedup of calculations by the SAM method and the BM$ method compared to the VΔC method. FIG. 10 is a diagram showing the results of an evaluation of the degree of speedup of vector-type arithmetic processing relative to scalar-type arithmetic processing. FIG. 10 illustrates an example of load balancing. FIG. 10 is a diagram showing the evaluation results of the degree of speedup of calculation processing by load distribution. FIG. 10 is a diagram illustrating an example of a measurement algorithm. FIG. 1 shows the relative speedup measured for the V.DELTA.C, SAM, and BM$ methods, and the memory hierarchy occupied depending on the problem size. FIG. 10 is a diagram illustrating an example of a ΔC generation circuit. FIG. 10 illustrates a first example of a hardware configuration for performing column swapping. FIG. 10 is a diagram illustrating an example of column interchange. FIG. 10 illustrates a second example of a hardware configuration for performing column swapping. FIG. 10 is a diagram illustrating a first modified example of the second example of the hardware configuration for performing column interchange. FIG. 10 is a diagram illustrating a second modification of the second example of the hardware configuration for performing column interchange. FIG. 10 illustrates a third example of a hardware configuration for performing column swapping. FIG. 10 is a diagram illustrating a modification of the third example of the hardware configuration for performing column interchange. FIG. 10 is a diagram illustrating another example of a ΔC generation circuit. FIG. 10 illustrates an example of a hardware configuration for performing processing on two replicas. FIG. 10 is a diagram illustrating another example of a ΔC generation circuit. FIG. 10 is a diagram illustrating an example of a replica processing circuit. FIG. 10 is a diagram illustrating an example of a replica processing circuit used in calculating a QAP using an asymmetric matrix. FIG. 1 is a diagram illustrating an example of hardware of a computer that is an example of a data processing device.

Hereinafter, embodiments of the invention will be described with reference to the drawings.
The data processing device of this embodiment searches for a solution to a quadratic assignment problem (QAP) or a quadratic semi-assignment problem (QSAP) by local search, as an example of an assignment problem. QAP, QSAP, and local search will be described below.

(QAP)
1 shows an example of how to calculate QAP. QAP is a problem of finding an allocation that minimizes the sum of the products of the flow rates between elements and the distances between the elements when allocating n elements (such as facilities) to n destinations.

FIG. 1 shows an example in which four facilities, each having identification numbers 1 to 4, are assigned to four locations (L ₁ to L ₄ ).
A flow matrix representing the amount of flow between n elements is expressed by the following equation (6).

The flow matrix (F) is a matrix with n rows and n columns. _{f i,j} is the flow amount in row i and column j, and represents the flow amount between elements with identification numbers = i, j. For example, the flow amount between a facility with identification number = 1 and a facility with identification number = 2 in Figure 1 can be expressed as f _1,2 .

The distance matrix representing the distance between n allocation destinations is expressed by the following equation (7).

The distance matrix (D) is a matrix with n rows and n columns. _{d k,l} is the distance in the kth row and lth column, and represents the distance between the assignees of identification numbers = k,l. For example, the distance between the assignee (L ₁ ) of identification number = 1 and the assignee (L ₂ ) of identification number = 2 in FIG. 1 can be expressed as d _1,2 .

The QAP is calculated by searching for an assignment that minimizes the above-mentioned equation (1). The assignment state of n elements to n assignment destinations is represented by an integer assignment vector φ or a binary state matrix X. φ is an element of the set Φ _n . Φ _n is a set of all permutations of the set N = {1, 2, 3, ..., n}. The binary variables x _i,j included in the binary state matrix X are expressed by the following equation (8).

1 shows an example in which a facility with identification number 1 is assigned to an assignee with identification number 2, a facility with identification number 2 is assigned to an assignee with identification number 3, a facility with identification number 3 is assigned to an assignee with identification number 4, and a facility with identification number 4 is assigned to an assignee with identification number 1. The integer assignment vector φ is φ(1) = 2, φ(2) = 3, φ(3) = 4, and φ(4) = 1. In the binary state matrix X, x _1,2 , x _2,3 , x _3,4 , and x _4,1 are each 1, and the other x _i,j are 0.

A QAP can have either one or both of the flow and distance matrices as symmetric matrices, or both asymmetric matrices. In this embodiment, we focus primarily on QAPs that use symmetric matrices (diagonal elements are 0 (biasless)). This is because such QAPs represent the majority of instances and simplify calculations. However, a QAP that uses a symmetric matrix can be directly converted to a QAP that uses an asymmetric matrix.

(QSAP)
QSAP is a modified version of QAP. In QSAP, the number of elements is not equal to the number of destinations. For example, multiple elements are allowed to be assigned to each destination. In QSAP, the distance matrix is expressed by the following equation (9):

The diagonal elements of the distance matrix (D) are non-zero to take into account routing within the assigned destination. QSAP also uses an additional matrix B, expressed as equation (10) below.

b _i,k represents a constant cost for assigning an element with identification number=i to an assignee with identification number=k.
The QSAP is calculated by searching for an allocation that minimizes the following equation (11) instead of equation (1), which is the evaluation function of the QAP.

The allocation state of n elements to m allocation destinations is expressed by an integer allocation vector ψ (ψε[1, m] ⁿ ) or a binary state matrix S. The binary variables s _i,j included in the binary state matrix S are expressed by the following equation (12).

(Local search)
In local search, a candidate solution is searched for within a nearby state that can be reached by changing the current state. A method for performing local search on a QAP is pairwise exchange between elements. In this method, two elements are selected and their assignment destinations are exchanged. The change in the value of the evaluation function (Equation (1)) due to the exchange of assignment destinations (represented by identification numbers = φ(a) and φ(b)) of two elements (represented by identification numbers = a and b) can be expressed by the following Equation (13).

As shown in equation (13), the change amount (ΔC _ex ) is generated by a multiply-and-accumulate loop.
When performing a local search on a QSAP by pairwise exchange, the amount of change in the value of the evaluation function (Equation (11)) can be expressed by the following Equation (14).

In equation (14), ΔB _ex indicates a change in a fixed cost related to the allocation of elements to each allocation destination, and is expressed by the following equation (15).

The state space in QSAP is not restricted to the permutations described above that represent the allocation state. An element can be relocated to a certain allocation destination regardless of whether another element is assigned to that allocation destination. The change in the value of the evaluation function when the element with identification number = a is reassigned from its current allocation destination (identification number = ψ(a)) to the allocation destination with identification number = l is expressed by the following equations (16) and (17).

In local search, a decision is made as to whether to accept a proposed allocation change that will cause a change in the value of the evaluation function by the amount of change calculated as above. The decision as to whether to accept a proposal is based on predefined criteria, such as a greedy method. When a greedy method is used, allocation changes that reduce the cost (value of the evaluation function (energy)) are accepted. When a greedy method is used, the proposal acceptance rate (PAR) is high at the start of the search, but tends to approach 0 later as the search gets stuck in a local minimum of the evaluation function and no further improvement is found.

Instead of the greedy method, in this embodiment, a stochastic local search method such as simulated annealing can be used. In the stochastic local search method, the acceptance probability (P acc ) of the Metropolis method, expressed by the following equation (18), can be used with a temperature parameter ( _T ) to add randomness to the change in allocation.

As T increases towards infinity, P _acc and PAR increase and all proposals are accepted regardless of the value of ΔC. As T is lowered and approaches 0, proposal acceptance becomes greedy and PAR tends to 0 as the search eventually gets stuck in a local minimum of the evaluation function.

The calculation of the QAP and ΔC of the QSAP takes up a large portion of the total processing time in the data processing device as it searches for a solution to the assignment problem.
The calculation of ΔC can be performed at high speed by using vector arithmetic operations as follows:

(Vectorization of ΔC calculation)
Hereinafter, this method may be referred to as the VΔC method. The calculation of ΔC _ex of QAP according to equation (13) can be easily vectorized using inner product calculations unless there is a condition that the calculation of i=a, b is skipped.

The data processing device of this embodiment eliminates this condition, performs the calculation of ΔC _ex as a product-sum operation loop for all n elements, and adds an additional product of the flow amount and the distance as a compensation term, as shown in the following equation (19), thereby performing compensation by eliminating the condition.

Equation (19) can be reformulated using the dot product, as in equation (22), using the flow rate difference vector ΔF (equation (20)) and the distance difference vector ΔD (equation (21)) when the allocation destinations of elements with identification numbers = a and b are swapped.

In equation (20), Δ ^b _a F represents the difference vector between row b and row a of the flow matrix, F _b,* represents row b of the flow matrix, and F _a,* represents row a of the flow matrix.
The elements of vector ΔF are arranged in the order of the original flow matrix, while the elements of vector ΔD are ordered to correspond to the current assignment state. This is mathematically shown by multiplying the transposed binary state matrix X (denoted as ^XT ) by (D _φ(a),* - D _φ(b),* ), as in equation (21). _{D φ(a),*} denotes row a of the distance matrix, and D _φ(b),* denotes row b of the distance matrix.

In software, calculations such as equation (21) are realized by rearranging the elements of vector ΔD in a time proportional to the number of elements=n.
In the case of QSAP, since a plurality of elements can be allocated to one allocation destination, a vector (ΔD) of size=n is generated from a distance matrix of m rows and m columns as shown in the following equation (23).

Using such a vector ΔD and the vector ΔF shown in equation (20), ΔC ^s _ex can be calculated by the following equation (24).

In the case of QSAP, since there is no restriction on multiple elements being assigned to the same destination, an additional compensation term such as the second line of equation (24) is added to ΔC.
Furthermore, the vectorized form of the calculation of ΔC ^s _rel by rearrangement corresponding to equation (17) can be expressed as the following equation (25).

Since the shift of the allocation destination of an element due to rearrangement is performed for one element, there is no need to calculate the vector ΔF and compensate by adding the product of the flow amount and the distance.
(Sorting the distance matrix)
The V.DELTA.C method can be implemented on processors such as multiple CPUs (Central Processing Units) and multiple GPUs (Graphics Processing Units) by using SIMD (Single Instruction/Multiple Data).

However, while the rearrangement of ΔD elements performed by the VΔC method can be performed in a time proportional to the number of elements (n), it is not very efficient because it is performed every time ΔC is calculated by swapping the allocation destinations of two elements. It can also incur a large amount of computational cost.

In order to minimize the computational cost of generating ΔD with an appropriate element order, the data processing device of this embodiment rearranges the columns of the distance matrix according to the current allocation state. Hereinafter, this method may be referred to as the SAM (State-Aligned D Matrix) method. Hereinafter, the distance matrix rearranged in this way is referred to as the state-aligned D matrix, and is denoted as D ^X when the assignment problem is QAP, and as D ^S when the assignment problem is QSAP.

FIG. 2 is a diagram showing an example of rearrangement of the distance matrix when calculating the QAP, and an example of a data processing device.
The data processing device 10 is, for example, a computer, and includes a storage unit 11 and a processing unit 12 .

The storage unit 11 is, for example, a volatile storage device that is an electronic circuit such as a dynamic random access memory (DRAM), or a non-volatile storage device that is an electronic circuit such as a hard disk drive (HDD) or flash memory. The storage unit 11 may also include an electronic circuit such as an SRAM (static random access memory) register.

The storage unit 11 stores, for example, a program that causes a computer to execute a local search or the like, as well as the flow matrix, distance matrix, state alignment D matrix (D ^X or D ^S ), and assignment state (represented by an integer assignment vector φ or a binary state matrix X) described above.

The processing unit 12 can be realized by a hardware processor such as a CPU, GPU, or DSP (Digital Signal Processor). The processing unit 12 may also be realized by an electronic circuit such as an ASIC (Application Specific Integrated Circuit) or FPGA. The processing unit 12 executes a program stored in the storage unit 11 to cause the data processing device 10 to perform local search processing. The processing unit 12 may also be a collection of multiple processors.

The processing unit 12 repeatedly switches the allocation destinations of the two elements, searching for an allocation state in which the value of the evaluation function shown in equation (1) or equation (11) is minimized, for example. The allocation state in which the smallest of the minimum values of the evaluation function is the optimal solution. Note that by changing the sign of the evaluation function shown in equation (1) or equation (11), the processing unit 12 can also search for an allocation state in which the value of the evaluation function is maximized (in this case, the maximum value is the optimal solution).

2 shows an example of rearrangement of a distance matrix when calculating a QAP. The allocation state at a certain time step = t and the state-ordered D matrix arranged corresponding to the allocation state are denoted as φ ^(t) and D ^X(t) . If an allocation change that swaps the allocation destinations of elements with identification numbers = 1 and 3 is proposed and accepted, the allocation state at time step = t+1 will be as shown by φ ^(t+1) in FIG. 2. In this case, D ^X(t+1) is generated as the state-ordered D matrix, in which the first and third columns of D ^X(t) are swapped to correspond to the accepted allocation state.

Note that if the proposed reassignment is not accepted, φ ^(t) and D ^X(t) are not updated.
In the case of QAP, the state alignment D matrix (D ^X ) can be expressed by multiplying the transpose matrix X ^T of the binary state matrix X representing the current allocation state by the original distance matrix (D), as shown in the following equation (26). Note that an example of a hardware configuration for rearranging the columns of the state alignment D matrix in hardware will be described later.

The vector ΔD ^X can be calculated by the following equation (27) without rearranging the elements of the above-mentioned vector ΔD, and ΔC _ex can be calculated by substituting ΔD ^X into the above-mentioned equation (22).

In the case of QSAP, the state alignment D matrix (D ^S ) can be expressed by multiplying the transpose matrix S ^T of the binary state matrix S representing the current allocation state by the original distance matrix (D) as shown in the following equation (28). ^{D S} is a matrix with m rows and n columns.

The vectors ΔD ^S resulting from the assignment exchange and rearrangement can be calculated using the following equations (29) and (30), respectively, without rearranging the elements of the above-mentioned ΔD, and by substituting these into the above-mentioned equations (24) and (25), ΔC ^s _ex and ΔC ^s _rel can be calculated.

In the VΔC method described above, rearranging the elements of vector ΔD requires operations proportional to the number of elements (n) each time ΔC is calculated. In contrast, in the above technique (SAM method), rearranging D ^X and D ^S to match the allocation state is only performed if the proposed allocation change is accepted. This reduces the number of operations required per iteration on average, improving computational efficiency, thereby shortening computation time and enabling the assignment problem to be calculated quickly.

(Boltzmann Machine Caching Method (Comparative Example))
The following method can also be considered as a method for vectorizing ΔC calculation. This method calculates and stores partial ΔC values (corresponding to the local field (h _i ) shown in the above-mentioned equation (3)) corresponding to bit inversions of the binary state matrix. Hereinafter, this method may be referred to as the Boltzmann machine caching method (BM$ method).

A cached local field is used for each bit of the binary state matrix (X or S). When a proposal for reassignment is accepted, an n-by-n matrix (hereinafter referred to as the cache matrix) using the local fields is updated in a time proportional to ^n² , but ΔC can be calculated in a time independent of n.

In the case of QAP, the cache matrix (H) is generated at the start of the search process according to equation (31) in a time proportional to n ³ .

Using the cached local field, the inner product ΔF·ΔD can be generated by equation (32). Then, by substituting the generated ΔF·ΔD into the above equation (22), ΔC _ex can be calculated.

If the reassignment proposal is accepted, the cache matrix is updated by using the Kronecker product according to equation (33), as shown in Figure 3 below.

FIG. 3 is a diagram illustrating an example of update calculation of the cache matrix.
The update is performed one row of the cache matrix at a time, skipping the processing of rows corresponding to zero elements of the vector ΔF.

In the case of QSAP, the cache matrix (H ^s ) is generated at the beginning of the search process according to equation (34).

Using the cached local field, the inner product ΔF·ΔD can be generated by equation (35). Then, by substituting the generated ΔF·ΔD into the above equation (24), ΔC ^s _ex can be calculated.

If the reassignment proposal is accepted, the cache matrix is updated using the Kronecker product according to equation (36).

Similarly, the cached local field can be used to generate the inner product ΔF·ΔD when rearrangement occurs using equation (37). Then, ΔC ^s _rel can be calculated by substituting the generated ΔF·ΔD into the above equation (25).

When a proposed relocation change is accepted, the cache matrix is updated according to equation (38).

In the BM$ method described above, a cache matrix is maintained, but as the problem size increases, the memory capacity required to maintain the cache matrix also increases. This makes it difficult to use relatively small memory that can be read quickly. The SAM method does not require the maintenance of such a cache matrix.

(Solver system design)
In order to compare and verify the performance of the above three methods (VΔC method, SAM method, and BM$ method), we designed the solver system shown below. Note that the following description will be given assuming that the processing unit 12 of the data processing device 10 includes a multi-core CPU with SIMD functionality and that the parallel tempering algorithm is implemented on the multi-core CPU.

Each dedicated core has a dedicated cache for holding search instance data (D ^x , D ^s , H, H ^s mentioned above).
FIG. 4 is a diagram illustrating an example of a solver system that performs parallel tempering.

The solver system has a search unit 20 and a parallel tempering controller 21. The search unit 20 has multiple cores 20a1-20am, each of which performs the aforementioned local search (SLS: Stochastic Local Search) on multiple replicas (corresponding to instances) in parallel. The number of replicas is M. Core 20a1 performs a local search on multiple replicas, including replica 20b1. Core 20am performs a local search on multiple replicas, including replica 20bM.

In the case of QAP, the aforementioned integer allocation vector φ, cache matrix (H) (in the case of the BM$ method), state alignment D matrix (D ^X ) (in the case of the SAM method), and value of the evaluation function (C) are stored for each replica.

In the case of QSAP, the aforementioned integer allocation vector ψ, cache matrix (H ^s ) (in the case of the BM$ method), state alignment D matrix (D ^s ) (in the case of the SAM method), and value of the evaluation function (C) are stored for each replica.

The flow matrix (F) shown in equation (6) and the matrix B shown in equation (10) used in the case of QSAP are common to all replicas 20b1 to 20bM.
The M replicas 20b1 to 20bM are set with different values of a temperature parameter (T). Each replica 20b1 to 20bM is set with one of temperatures that increases stepwise from T _min to T _max . Hereinafter, such stepwise increasing temperature parameter values may be referred to as a temperature ladder.

The parallel tempering controller 21 exchanges the values of Tk and _{Tk+1 between replicas in which adjacent temperature parameter values (Tk, Tk+1 )} are set, based on the evaluation function values ( _Ck , _{Ck +1} ) and _Tk , Tk _{+ 1} , with the exchange allowance probability SAP expressed by the following equation (39).

The search unit 20 and the parallel tempering controller 21 described above are realized by the processing unit 12 and the storage unit 11 shown in FIG.
The flow matrix, state alignment D matrix, and cache matrix are stored in the storage unit 11 in the single-precision floating-point format, for example, so that a fused multiply-and-accumulate instruction can be used when performing inner product calculations when executing the VΔC method or the SAM method, or when updating the cache matrix.

Figure 5 shows an example of an algorithm for finding a solution to QAP using local search with parallel tempering. In the example algorithm shown in Figure 5, the number of replicas for which local search is performed in parallel is 32. In Figure 5, "(1)" and "(31)" on lines 16 and 18 represent the aforementioned equations (1) and (31).

The algorithm begins by initializing the temperature parameter values and replica states. The solver system then enters an optimization loop, performing one local search iteration in parallel across all replicas.

After I local searches, a replica exchange process is performed and the optimization loop is restarted until a preset Best-Known-Solution (BKS) cost (C _BKS ) is found or a timeout limit is reached.

In the local search of the QAP, a loop is executed for a preset number of iterations (I), two elements (facilities in the example of FIG. 5) are selected, and the value of ΔC _ex when the allocation destinations of the two elements are swapped is calculated using the selected method (VΔC method, BM$ method, or SAM method).

Next, P _acc in equation (18) is calculated using the calculated ΔC _ex and the current temperature parameter value of each replica. A random value between 0 and 1 is then generated and compared with P _acc to perform a Bernoulli trial to determine whether the proposed swap of assignments between the two selected elements is acceptable. If the proposal is accepted, the state (assignment state) is updated.

Figure 6 is a flowchart showing the overall processing flow of local search using parallel tempering. Figure 6 shows the overall processing flow of local search using parallel tempering for the algorithm shown in Figure 5.

First, an initialization loop (steps S10 to S13) is performed while incrementing i by 1 from i=0 to i<M-1, where M is the number of replicas (32 in the example of FIG. 5).
In the initialization loop, the initial temperature (temperature ladder) is set for each replica (T[i]←T ⁰ [i]) (step S11), and then the replica is initialized (see FIG. 7) (step S12).

Next, an optimization loop (steps S14 to S16) is performed while increasing i by one from i=0 to i<M-1.
In the optimization loop, a replica search process is performed (step S15), in which a local search is performed I times for each replica using the set temperature parameter value.

Then, a replica exchange process (step S17) is performed, and it is determined whether or not to terminate the search (step S18). For example, as described above, if a preset BKS cost (C _BKS ) is found or the timeout limit is reached, it is determined that the search is terminated, and the search ends. If it is determined not to terminate the search, the optimization loop is restarted from step S14.

When the search is completed, the data processing device 10 may output the search results (e.g., minimum values C _min and φ _min described below). The search results may be displayed on a display device connected to the data processing device 10, or may be transmitted to an external device.

Note that while Figure 6 shows an example in which processing is performed sequentially for each replica, it is possible to execute steps S11, S12, and S15 in parallel for 32 replicas using a processor with 32 cores, for example.

Next, the replica initialization process in step S12 and the replica search process in step S15 will be described with reference to a flowchart, taking the case where the SAM method is executed as an example.
FIG. 7 is a flowchart showing an example of the flow of the replica initialization process in the case of a QAP.

First, the integer allocation vector φ is initialized (step S20). In step S20, the allocation destination of elements (location of facilities) is determined randomly, for example. Then, the initial cost (C) is calculated using equation (1) (step S21).

Then, the minimum values (C _min and φ _min ) are initialized by the initialized φ and C (step S22). Also, the state alignment D matrix (D ^X ) is initialized to correspond to the initial value of φ (step S23). Then, the process returns to the flowchart shown in FIG.

FIG. 8 is a flowchart showing an example of the flow of replica search processing in the case of QAP.
In the replica search process, an iteration loop (steps S30 to S40) is performed while increasing i by one from i=1 to i<I.

First, two elements (represented by identification numbers a and b) that are candidates for replacement of allocation destinations are selected (step S31).
Then, the value of ΔC _ex (a, b) when the allocation destinations of the two elements are swapped is calculated using equation (22) (step S32).

Next, P _acc in equation (18) is calculated using the calculated ΔC _ex and the current temperature parameter value of each replica (step S33), and it is determined whether P _acc is greater than a random number rand() ranging from 0 to 1 (step S34).

If it is determined that P _acc > rand(), columns a and b of D ^X are swapped (D ^X _{*, a} , D ^X _{*, b} ←D ^X _{*, b} , D ^X _{*, a} ) (step S35). Furthermore, the allocation status is updated by swapping the allocation destinations (φ(a), φ(b) ←φ(b), φ(a)) (step S36), and the value of the evaluation function (C) is updated (C ← C + ΔC _ex (a, b)) (step S37).

Thereafter, it is determined whether C<C _min (step S38), and if it is determined that C<C _min , the minimum value is updated (C _min ←C, φ _min ←φ) (step S39).

After the process of step S39, or if it is determined in the process of step S34 that P _acc > rand() is not true, or if it is determined in the process of step S38 that C < C _min is not true, the process from step S31 is repeated until i = I. If i = I, the process returns to the process of the flowchart shown in FIG.

The order of the processes shown in FIGS. 6 to 8 is an example, and the order of the processes may be changed as appropriate.
FIG. 9 is a diagram showing an example of an algorithm for searching for a solution to QSAP by local search using parallel tempering.

The local search function for computing the QSAP shown in Figure 9 includes a rearrangement loop that determines whether to accept or update a rearrangement proposal, before an exchange loop that determines or updates whether to swap the allocation destinations of two elements.

Next, using the SAM method as an example, we will explain the replica initialization process in step S12 of Figure 6 and the replica search process in step S15 when calculating QSAP using a flowchart.

FIG. 10 is a flowchart showing an example of the flow of a replica initialization process in the case of QSAP.
First, the integer allocation vector ψ is initialized (step S50). In the processing of step S50, the allocation destinations of elements (location destinations of facilities) are determined, for example, randomly. Then, the initial cost (C) is calculated using equation (11) (step S51).

Then, the minimum values (C _min and ψ _min ) are initialized by the initialized ψ and C (step S52). Also, the state alignment D matrix (D ^S ) is initialized by equation (28) so as to correspond to the initial value of ψ (step S53). Then, the process returns to the flowchart shown in FIG.

FIG. 11 is a flowchart showing an example of the flow of replica search processing in the case of QSAP.
In the replica search process in the case of QSAP, an iteration loop (steps S60 to S78) is also performed while increasing i by one from i=1 to i<I.

First, an element with identification number = a and an allocation destination of identification number = l are selected (step S61).
Then, the value of ΔC ^s _rel when the element with identification number=a is assigned to the assignment destination with identification number=l is calculated using equation (25) (step S62).

Next, P _acc in equation (18) is calculated using the calculated ΔC ^s _rel and the current temperature parameter value of each replica (step S63), and it is determined whether P _acc is greater than a random number rand() ranging from 0 to 1 (step S64).

If it is determined that P _acc > rand(), the values of column a of D ^s are updated with the values of column l of the distance matrix D (step S65). Furthermore, the allocation state is updated (φ(a)←l) (step S66), and the value of the evaluation function (C) is updated (C←C+ΔC ^s _rel ) (step S67).

Thereafter, it is determined whether C<C _min (step S68), and if it is determined that C<C _min , the minimum value is updated (C _min ←C, ψ _min ←ψ) (step S69).

After the process of step S69, or if it is determined in the process of step S64 that P _acc > rand( ) is not true, or if it is determined in the process of step S68 that C < C _min is not true, the process of step S70 is performed.

In the process of step S70, two elements (represented by identification numbers a and b) are selected as candidates for replacement of allocation destinations.
Then, the value of ΔC ^s _ex when the allocation destinations of the two elements are swapped is calculated using equation (24) (step S71).

Next, it is determined whether ΔC ^sex _< 0 (step S72). Note that a predetermined value (fixed value) may be used instead of 0.
If it is determined that ΔC _ex ^s < 0, columns a and b of D ^S are swapped (D ^S _*,a , D ^S _*,b ←D ^S _*,b , D ^S _*,a ) (step S73). Furthermore, the allocation state is updated by swapping the allocation destinations (ψ(a), ψ(b) ←ψ(b), ψ(a)) (step S74), and the value of the evaluation function (C) is updated (C ←C + ΔC ^s _ex ) (step S75).

Thereafter, it is determined whether C<C _min is true (step S76), and if it is determined that C<C _min is true, the minimum value is updated (C _min <-C, ψ _min <-ψ) (step S77).

After the process of step S77, or if it is determined in the process of step S72 _that ΔC ^sex < 0 is not satisfied, or if it is determined in the process of step S76 that C < C _min is not satisfied, the process from step S61 onwards is repeated until i = I. If i = I, the process returns to the process of the flowchart shown in FIG.

The order of the processes shown in FIGS. 10 and 11 is an example, and the order of the processes may be changed as appropriate.
(Evaluation results of calculation speed)
First, we will show the results of evaluating the calculation speed of three methods using scalar arithmetic processing: the VΔC method, the SAM method, and the BM$ method. The calculation targets 10 QAP instances, and calculations were performed using the algorithm shown in Figure 5. The calculations were performed using two solver systems such as those shown in Figure 4, with 64 replicas in parallel tempering.

Each instance was run sequentially 100 times with independent random seeds, and the time to reach BKS (Time-to-Solution (TtS)) was recorded.
FIG. 12 shows the evaluation results of the degree of speedup of calculations by the SAM method and the BM$ method compared to the VΔC method. The horizontal axis represents 10 instances of QAP and the geometric mean (denoted as "GEOMEAN" in FIG. 12), and the vertical axis represents the degree of speedup of the SAM method and the BM$ method compared to the VΔC method. The degree of speedup is represented by the ratio of TtS by the SAM method and the BM$ method to TtS by the VΔC method.

When the SAM method is used, the calculation speed is improved compared to when the VΔC method is used, despite the additional calculation cost of updating ^DX .
When the BM$ method is used, the calculation speed is improved by 2.57 times in geometric mean compared to the VΔC method.

(Improvement of calculation speed in the case of SIMD)
Next, we will show the results of evaluating the degree of improvement in calculation speed of vector-type arithmetic processing compared to scalar-type arithmetic processing. To perform the vector-type arithmetic processing, SIMD built-in functions of AVX (Advanced Vector eXtensions) 2 were used.

Using the same procedure as above, the same 10 QAP instances were run with vector-based computation.
13 is a diagram showing the evaluation results of the degree of speedup of vector-type arithmetic processing relative to scalar-type arithmetic processing. The horizontal axis represents 10 instances of QAP and the geometric mean (denoted as "GEOMEAN" in FIG. 13), and the vertical axis represents the degree of speedup of vector-type arithmetic processing relative to scalar-type arithmetic processing. The degree of speedup is represented by the ratio of TtS by vector-type arithmetic processing to TtS by scalar-type arithmetic processing for each of the VΔC method, the SAM method, and the BM$ method.

On average, the VΔC method is nearly twice as fast as scalar-type arithmetic processing in vector-type arithmetic processing, and the BM$ and SAM methods are more than three times faster. As mentioned above, the VΔC method offers the least benefit from SIMD built-in functions, since the rearrangement of the ΔD elements performed in the VΔC method takes up a significant portion of the execution time.

(Dynamic Load Balancing)
In parallel tempering, the PAR differs between the replica with the lowest temperature setting and the replica with the highest temperature setting. The larger the value of T (the higher the temperature), the higher the PAR. This can cause large execution time gaps between the threads processing the replicas in the SAM and BM$ update processes.

To mitigate this, the data processing device 10 tracks the time per iteration at each temperature within the SLS function and uses that time to scale the number of iterations performed at each temperature.

Figure 14 shows an example of load balancing. Figure 14 shows the calculation time of ΔC before exchanging temperature parameters and the update processing time for 32 replicas set with a temperature ladder of T1 to T32 (T1<T2<...<T32).

Without load balancing, replicas with higher temperature parameter values take longer to update, while replicas with lower temperature parameter values take shorter update times. As a result, replicas with lower temperature parameter values experience longer idle times.

When load balancing is performed, for example, the number of iterations of the ΔC calculation in other replicas is scaled based on the execution time (thread execution time) of the replica to which T1 is set.

This allows the thread execution times to be roughly equal among the replicas, as shown in FIG. 14, thereby reducing idle time.
In order to quantitatively demonstrate the effectiveness of such load balancing, the degree of speedup of calculation processing due to load balancing was evaluated using the above-mentioned 10 instances as calculation targets, and the results are shown below.

Figure 15 shows the evaluation results of the degree of speedup of calculation processing due to load balancing. The horizontal axis represents 10 instances of QAP and the geometric mean "GEOMEAN", and the vertical axis represents the degree of speedup when load balancing is performed compared to calculation processing without load balancing. The degree of speedup is represented by the ratio of TtS when the SAM method and the BM$ method are performed with load balancing to TtS when the SAM method and the BM$ method are performed without load balancing. Note that the same set of temperature parameter values (for example, T1 to T32 in Figure 14) was used whether load balancing was performed or not.

We found that load balancing achieved an average speedup of 1.86 times for the SAM method and 1.38 times for the BM$ method. The difference in the effectiveness of load balancing between instance families may be due to different PAR gaps between extreme values of the temperature parameter, depending on the characteristics of the instance families.

Table 1 summarizes the evaluation results of the computational speedup provided by each of the above functions.

Two degrees of speed-up are shown (Incremental Speed-Up and Cumulative Speed-Up). Note that the degree of speed-up is calculated using the geometric mean of the TtS for the 10 instances mentioned above, with the TtS obtained by scalar-type calculation processing using the VΔC method set as the standard (1.00).

(Benchmark results)
Regarding the V.DELTA.C method, the SAM method, and the BM$ method, benchmark scores were measured using a data processing device 10 with a predetermined hardware configuration including a 64-core CPU.

To measure the benchmark scores of QAP solvers that implement the above three methods, instances from the QAP Library Reference (see Non-Patent Document 5) were used. Additionally, instances proposed in Non-Patent Documents 6 and 7 were used. In this embodiment, we also attempted to find a new BKS state for the large-size QAP (with no known optimal solution) shown in Non-Patent Document 7. Furthermore, to measure the benchmark scores of QSAP solvers, a set of instances introduced in Non-Patent Document 8 was used.

(QAP benchmark)

Table 2 shows the TtS for each instance of the QAP using solvers that implement three vector-based computation methods (the VΔC method, the SAM method with load balancing, and the BM$ method). Table 2 also shows the degree of speedup (calculated using the geometric mean of TtS), with the TtS for vector-based computation using the VΔC method set as the baseline (1.00). Table 2 also shows the TtS and degree of speedup for two publicly available solvers, ParEOTS (see Non-Patent Document 9) and PBM (Permutational Boltzmann Machine) (see Non-Patent Document 1). These two publicly available solvers are capable of solving some difficult instances with a 100% success rate within a 5-minute timeout window.

The TtS values shown in Table 2 are the average TtS values (expressed as values with a 99% confidence interval) measured over 100 independent consecutive runs of local searches using the VΔC, SAM, and BM$ methods. The timeout limit for each run is 5 minutes. The TtS values for ParEOTS and PBM are from their respective publications.

As shown in Table 2, the BM$ method performed best across all instances among the five solvers, achieving a speedup of 2x compared to PBM and over 300x compared to ParEOTS. The SAM and BM$ methods achieved average speedups of 1.92x and 3.22x compared to the VΔC method, respectively. Note that the differences in the results in Table 2 compared to Table 1 are due to the instances used.

However, whether the SAM or BM$ method is better depends on the problem scale and PAR, and there are cases where the SAM method is superior to the BM$ method. For example, as will be described later (see Figure 17), as the PAR increases, the SAM method tends to be superior to the BM$ method. Furthermore, these evaluation results are the result of an implementation on a CPU, and are the results of a system with a relatively large memory. When implemented on a dedicated circuit, the SAM method, which does not need to store local fields, may be superior to the BM$ method in some cases.

(QSAP Benchmark)
Compared to QAP, there has been little previous disclosure on QSAP solvers. The main reference is the Parallel Memetic Iterative Tabu Search (PMITS) algorithm implemented on a 20-core CPU (see Non-Patent Document 8). PMITS reaches the best solution with a 1-hour timeout limit using 50% of the population, which is the same as the stopping criterion. This method is a reasonable way to predict convergence in memetic algorithms with cooperative replicas. However, in parallel tempering, where replicas are constantly flowing along the temperature ladder, such a criterion is not valid.

Therefore, in this embodiment, TtS was measured in the same manner as for QAP, and the TtS value was measured as a reference for PMITS. The results are shown in Table 3 below. No correlation between the parallel tempering solver and PMITS is shown.

As shown in Table 3, when compared to QSAP, the performance of the SAM and BM$ methods is almost 2x worse than the VΔC method across QSAP instances. This is likely due to higher PAR.

QAP Scaling on Extended Taillard Instances
Some QAP instances introduced in Non-Patent Document 7 have unknown optimal solutions and have a maximum size of n = 729. Due to their size and difficulty, they are rarely used for benchmarking, but the data processing device 10 executed these instances using the VΔC method, the SAM method, and the BM$ method with a predetermined time limit to find better solutions.

Previous attempts to improve the BKS values of these instances have resulted in little or no improvement, even after running for minutes to hours (see, for example, Non-Patent Document 10).

In the experiment, each instance was executed 20 times, completing in 10 seconds for n=125 and n=175, and 30 seconds for n=343 and n=729 instances. Tables 4 and 5 show the average cost along with the best cost (the value of the evaluation function in Equation (1)) when applying the VΔC, SAM, and BM$ methods for each instance.

Despite the shorter execution time, we were able to improve BKS for all but four instances using the SAM and BM$ methods, and found that the average cost was also lower than the previous BKS.

(Scaling analysis)
Based on benchmark results and a qualitative comparison of the VΔC and SAM methods, the SAM method appears to be more efficient because the VΔC method sequentially reorders each element of ΔD at each iteration, whereas the SAM method reorders D ^X and D ^S to match the allocation state only if the proposed reassignment is accepted.

On the other hand, the relative performance of SAM and BM$ methods depends primarily on PAR, which can vary significantly depending on the search algorithm used. Furthermore, even within the same search algorithm, PAR can change between runs using algorithms such as simulated annealing, or between instances searched simultaneously, such as in parallel tempering.

Below, we show the results of comparing execution times according to problem size and PAR value.
FIG. 16 is a diagram illustrating an example of a measurement algorithm.
In the measurement algorithm, first, a flow matrix (F) and a distance matrix (D) with random values between 1 and 100 are generated. Then, the allocation state (φ) is initialized. Next, a predetermined number of iterations (I _limit ) are executed, and the execution time of one loop is measured. At this time, a proposal with a desired PAR value is randomly accepted.

The problem sizes are divided into three groups based on the memory hierarchy used to store replica data for the SAM and BM$ methods. For problem sizes of n<256 and 256<n<1024, I _limit is set to 100M and 10M, respectively. For problem sizes of n>1024, I _limit is set to 1M.

At each iteration, the value of ΔC _ex is calculated and a random number generator is used to determine whether to accept the proposal with the desired PAR.
To measure performance while the CPU pipeline, cache, and memory were fully loaded, 64 loop instances were run in parallel (one per core). This process was repeated using 10 different random number seeds for each data point, and the average execution time was measured.

Because the density of the flow matrix affects the update function of the cache matrix (H) in the BM$ method, two separate sets of simulations were run: one with a fully dense flow matrix and one with a sparse flow matrix with only 10% non-zero values. Each measurement was calculated for 10 different problem sizes in the range n = [100, 5000] and 29 PARs in the range PAR = [0.0001, 0.1], for a total of 290 parameter combinations.

Figure 17 shows the relative speedup measured for the VΔC, SAM, and BM$ methods, and the memory hierarchy occupied depending on the problem size. In the diagram showing the speedup, the horizontal axis represents PAR [%], and the vertical axis represents the speedup (Speed-Up). Figure 17 also shows the speedup of the BM$ method relative to the SAM and VΔC methods, using flow matrices of different densities. Furthermore, Figure 17 shows the speedup of the SAM method relative to the VΔC method. In the example of Figure 17, the memory hierarchy is arranged in order of decreasing storage capacity: L2 cache, L3 cache, and DRAM.

Note that the maximum PAR values for a given flow matrix density in a non-load-balanced parallel tempering simulation for 10 instances of QAP were as shown in Table 6 below.

The QSAP simulation was omitted because it is similar to the QAP simulation, and it is believed that there is no significant difference between the results of both simulations.
The relative speedups for different problem sizes can be divided into three groups based on which layer of the memory hierarchy contains the majority of the data used for searching: Each core has its own L2 cache (256 kB in this example), and a group of four cores shares an L3 cache (16 MB in this example).

As shown in FIG. 17, in the case of the SAM and BM$ methods, the D ^X and cache matrix (H) of each loop thread can fit into the L2 cache of the core with a maximum size of n=256, and problems up to n=1024 can fit into the L3 cache.

As shown in Figure 17, as the problem size increases and a memory hierarchy with larger storage capacity is used, the relative performance (degree of speedup) between the BM$ method and the other two methods decreases. As the data used for the search moves to a higher memory hierarchy, the PAR value required for the BM$ method to maintain superior performance to the SAM and VΔC methods decreases significantly.

For solvers that implement parallel tempering, the maximum PAR across replicas does not necessarily decrease with problem size across instance families, as shown in Table 6. This indicates that for smaller QAP instances, the relationship between PAR and problem size to maintain the same relative speeds as shown in Table 2 does not necessarily depend on the relationship shown in Figure 17.

The speedup of the SAM method over the VΔC method can be divided into two groups based on the problem size. When the problem size is n≦800 and most of the D ^X for each loop instance fits in the cache, the SAM method is clearly more advantageous than the VΔC method in the range of PAR<10%. Moving the storage destination of the search data from the L2 cache to the L3 cache does not have a significant impact on the relative performance between the two methods.

(Hardware example)
An example of hardware for implementing the SAM method will be described below.
In the following, for the sake of simplicity, it is assumed that both the flow matrix and the distance matrix are symmetric matrices (diagonal elements are 0 (biasless)). As mentioned above, this type of QAP is the majority of instances, and it simplifies calculations. A QAP using a symmetric matrix can be directly converted into a QAP using an asymmetric matrix.

The evaluation function for QAP using only symmetric matrices can be expressed as follows:

The only difference between equation (1) and equation (40) is that "j=1" in equation (1) becomes "j=i".
When the evaluation function expressed by equation (40) is used, ΔC _ex used in calculating QAP can be expressed by the following equations (41) and (42) instead of equations (19) and (22).

FIG. 18 is a diagram illustrating an example of a ΔC generating circuit.
The ΔC generation circuit 30 includes a flow matrix memory 31 a, a state-ordered D matrix memory 31 b, difference calculation circuits 32 a and 32 b, multiplexers 33 a and 33 b, an inner product calculation circuit 34, a register 35, a multiplication circuit 36, and an addition circuit 37. These are, for example, circuits or memories included in the storage unit 11 or the processing unit 12 shown in FIG.

The flow matrix memory 31a stores the flow matrix (F).
The state alignment D matrix memory 31b stores the state alignment D matrix (D ^X ) which is an updated distance matrix.

In the example of Figure 18, the flow matrix memory 31a and the state-aligned D matrix memory 31b are dual-port memories with two ports. Such dual-port memories can be used in some FPGAs.

The difference calculation circuit 32a calculates F _b,* - _{F a,*,} which is Δ ^b _a F in equation (42).
The difference calculation circuit 32b calculates D ^{X φ(a)} ,* - D ^X _φ(b) ,*, which is Δ _φ(a) ^φ _(b) D ^X in equation (42).

The multiplexer 33a selects and outputs f _a,b from F _a,* .
The multiplexer 33b selects and outputs d _φ(a),b from D ^X _φ(a),* .
The dot product calculation circuit 34 calculates the dot product of Δ ^b _a F and Δ ^φ(a) _φ(b) D ^X. The dot product calculation circuit 34 can be realized by, for example, a plurality of multipliers connected in parallel.

The register 35 holds the coefficient "2" included in the equation (42).
The multiplication circuit 36 calculates 2f _a,b d _φ(a),b .
The adder circuit 37 calculates and outputs ΔC _ex by adding 2f _a,b d _φ(a),b to the calculation result of the inner product.

The calculation of ΔC _ex using such hardware is controlled by, for example, a control circuit (not shown) included in the processing unit 12 of FIG. 2 executing a program.
As described above, in the SAM method, when a swap (allocation change) of the allocation destination of the element that generates ΔC _ex is proposed and accepted, the columns of the state alignment D matrix are swapped to correspond to the accepted assigned state.

FIG. 19 illustrates a first example of a hardware configuration for performing column interchange.
19, the columns of the state order D matrix stored in the state order D matrix memory 31b can be interchanged using, for example, multiplexers 40a and 40b and a switch 41. These circuits may also be included in the processing unit 12 in FIG.

The multiplexer 40a sequentially selects and outputs the value of the first column of each row of the state alignment D matrix read out from the state alignment D matrix memory 31b for replacement.
The multiplexer 40b sequentially selects and outputs the value of the second column of each row of the state alignment D matrix read from the state alignment D matrix memory 31b.

Switch 41 swaps the memory locations in state-aligned D-matrix memory 31b so that the value output from multiplexer 40a is stored in the location where the value output from multiplexer 40b was stored. Switch 41 also swaps the memory locations in state-aligned D-matrix memory 31b so that the value output from multiplexer 40b is stored in the location where the value output from multiplexer 40a was stored.

Such column switching using hardware is controlled by, for example, a control circuit (not shown) included in the processing unit 12 of FIG. 2 executing a program.
20 is a diagram showing an example of column interchange, in which the first and third columns are interchanged.

First, d _1,3 and d _1,1 are selected by multiplexers 40a and 40b, and their storage locations are swapped by switch 41. Next, d _2,3 and d _2,1 are selected by multiplexers 40a and 40b, and their storage locations are swapped by switch 41. The same process is repeated a total of n cycles to complete the column swapping.

Fig. 21 is a diagram showing a second example of a hardware configuration for performing column interchange, in which the same elements as those shown in Fig. 19 are denoted by the same reference numerals.
21 , the state order D matrix memory 31b stores not only the state order D matrix but also the transposed matrix ((D ^X ) ^T ) of the state order D matrix. The columns of the state order D matrix can be swapped using the elements of the read transposed matrix, for example, by the aforementioned switch 41 and shift registers 45a and 45b. These circuits may also be included in the processing unit 12.

In the circuit configuration of FIG. 21, two rows of the transposed matrix corresponding to the two columns of the state reordering D matrix to be transposed are read out from the state reordering D matrix memory 31b.
The shift register 45a holds n values in the first row of the two rows of the transposed matrix read from the state-ordered D-matrix memory 31b, and shifts them one cycle at a time to output the values one by one.

Shift register 45b holds the n values in the second row of the two rows of the transposed matrix read from state-aligned D-matrix memory 31b, shifts them one cycle at a time, and outputs each value one by one.

Switch 41 switches the storage location in state-ordered D-matrix memory 31b so that the value output from shift register 45a is stored in the location where the value output from shift register 45b was stored. Switch 41 also switches the storage location in state-ordered D-matrix memory 31b so that the value output from shift register 45b is stored in the location where the value output from shift register 45a was stored.

This hardware configuration also allows columns to be swapped every n cycles. Furthermore, the hardware configuration of Figure 21 may have improved wiring efficiency compared to the hardware configuration of Figure 19. It also eliminates the need for multiplexers 40a and 40b as shown in Figure 19, which may result in a relatively large circuit size.

Figure 22 shows a first variant of the second example of a hardware configuration for column swapping. In Figure 22, elements that are the same as those shown in Figure 21 are designated by the same reference numerals.

The hardware configuration in Fig. 22 differs from the hardware configuration in Fig. 21 in that the transpose matrix ((D ^X ) ^T ) of the state ordered D matrix is stored in a transpose matrix memory 46 that is separate from the state ordered D matrix memory 31b. The above two rows of the transpose matrix are read from the transpose matrix memory 46. The other configurations are the same as those in Fig. 21.

Figure 23 shows a second variation of the second example of a hardware configuration for column swapping. In Figure 23, elements that are the same as those shown in Figure 22 are designated by the same reference numerals.

In the hardware configuration of Fig. 23, in addition to the flow matrix, the transposed matrix ((D ^X ) ^T ) of the state order D matrix is also stored in the flow matrix memory 31a shown in Fig. 18. The columns of the state order D matrix can be swapped by the switch 41 and shift registers 45a and 45b, as described above, using the elements of the transposed matrix read from the flow matrix memory 31a.

Fig. 24 is a diagram showing a third example of a hardware configuration for performing column interchange, in which the same elements as those shown in Fig. 19 are denoted by the same reference numerals.
In the hardware configuration of FIG. 24, the state alignment D matrix is stored in two state alignment D matrix memories 31b1 and 31b2.

Using the elements of each row of the state alignment D matrix read from state alignment D matrix memory 31b1, columns of the state alignment D matrix stored in state alignment D matrix memory 31b2 are swapped in the same way as when the hardware configuration of Figure 19 is used.

The state-ordered D matrix after the column exchange is copied to the state-ordered D matrix memory 31b1.
With this hardware configuration, columns can be swapped in n cycles. In addition, the hardware configuration of Fig. 24 may have improved wiring efficiency compared to the hardware configuration of Fig. 19.

25 is a diagram showing a modification of the third example of the hardware configuration for performing column interchange, in which the same elements as those shown in FIG. 19 are denoted by the same reference numerals.
In the hardware configuration of Fig. 25, in addition to the flow matrix, the transposed matrix ((D ^X ) ^T ) of the state order D matrix is also stored in the flow matrix memory 31a shown in Fig. 18. The columns of the state order D matrix can be swapped by the multiplexers 40a and 40b and the switch 41, as described above, using the elements of the transposed matrix read from the flow matrix memory 31a.

Fig. 26 is a diagram showing another example of a ΔC generating circuit, in which the same elements as those shown in Fig. 18 are denoted by the same reference numerals.
26, single-port memories are used as the flow matrix memory 51a and the state-ordered D matrix memory 51b. In this case, a register 52a is used that holds F _a,* , outputs F _a, * when F _b,* is read from the flow matrix memory 51a, and supplies it to the difference calculation circuit 32a and the multiplexer 33a. Also, a register 52b is used that holds D ^x _φ(a),* , outputs D ^x φ(a),* when D ^x φ _{( b),*} is read from the state-ordered D matrix memory 51b, and supplies it to the difference calculation circuit 32b and the multiplexer 33b.

The other configurations are the same as those in FIG.
Incidentally, when the hardware configuration for column exchange shown in FIGS. 21 to 24 is applied to the ΔC generation circuits 30 and 50 shown in FIG. 18 or 26, it becomes possible to simultaneously execute processing for two replicas.

Figure 27 shows an example of a hardware configuration for processing two replicas. Figure 27 shows a configuration that combines the ΔC generation circuit 50 shown in Figure 26 with the hardware configuration for column swapping shown in Figure 22.

In the example of FIG. 27, the state alignment D matrix memory 51b stores state alignment D matrices (D ^x _R1 , D ^x _R2 ) of two replicas (R1, R2).
One replica accepts the assignment change and performs column swapping for that replica using the write port of the state-aligned D-matrix memory 51 b, while the other replica performs a process of calculating ΔC _ex based on the elements of the state-aligned D-matrix read from the read port of the state-aligned D-matrix memory 51 b.

Note that if neither replica is currently performing an update, processing for one replica will be stalled while ΔC _ex is being calculated for the other replica.

Fig. 28 is a diagram showing another example of a ΔC generating circuit, in which the same elements as those shown in Fig. 18 are denoted by the same reference numerals.
When only the problem of symmetric permutation is the object of calculation and ΔC _ex is expressed by the following equation (43), a ΔC generating circuit 60 as shown in FIG. 28 can also be used.

In the ΔC generating circuit 60, the multiplexers 33a and 33b shown in FIG. 18 are not necessary, and instead a selector 61 is provided.
For example, the difference calculation circuit 32b has n subtractors 32bi that calculate d _φ(a),i −d _φ(b),i , and the selector 61 has n two-input, one-output multiplexers 61ai that select and output either the output of the subtractor 32bi or 0.

Multiplexer 61ai outputs 0 when i = a, b in equation (43). In such cases, the method of outputting 0 may be other than using a multiplexer.

(Replica processing circuit)
The circuit that performs processing for each replica can be realized by combining any of the ΔC generation circuits 40, 50, and 60 described above with any of the several hardware configurations described above that perform column interchange of the state alignment D matrix.

Figure 29 is a diagram showing an example of a replica processing circuit. Figure 29 shows an example of a replica processing circuit 70 that combines the ΔC generation circuit 50 shown in Figure 26 with the hardware configuration for column swapping shown in Figure 19. In Figure 29, elements that are the same as those shown in Figures 19 and 26 are assigned the same reference numerals.

In the example of FIG. 29, the multiplexer 33b also has the function of the multiplexer 40a shown in FIG.
These circuits may be included in the memory unit 11 or the processing unit 12 in FIG.

Note that the configuration for storing the integer allocation vector φ representing the allocation state, and the configuration for determining whether or not to accept a proposal for an allocation change that will result in the calculated ΔC _ex are not shown in the figure.

(Extension to QAP when using asymmetric matrices)
In the calculation of QAP using an asymmetric matrix (diagonal elements are non-zero), ΔC _asym , which corresponds to ΔC _ex when a symmetric matrix is used, is expressed by the following equation (44).

Such a replica processing circuit for calculating ΔC _asym can be realized using a replica processing circuit for calculating ΔC _ex .
30 is a diagram showing an example of a replica processing circuit used in calculating QAP using an asymmetric matrix. In FIG. 30, elements that are the same as those shown in FIG. 29 are assigned the same reference numerals.

Replica processing circuit 80 that calculates ΔC _asym has two replica processing circuits 70a1 and 70a2 that calculate ΔC _ex shown in Fig. 29. However, the flow matrix memory 51a of one of the replica processing circuits, 70a2, stores the transposed matrix (F ^T ) of the flow matrix, and the state order D matrix memory 51b of replica processing circuit 70a2 stores the transposed matrix ((D ^X ) ^T ) of the state order D matrix.

The replica processing circuit 80 further includes memories 81a and 81b, registers 82a and 82b, difference calculation circuits 83a and 83b, a multiplication circuit 84, a compensation term calculation circuit 85, and an addition circuit 86.

The memory 81a stores the diagonal elements of the flow matrix ( _Fd ), and the memory 81b stores the diagonal elements of the distance matrix ( _Dd ). When an asymmetric matrix is used, unlike when a symmetric matrix is used, these diagonal elements may contain non-zero values.

The register 82a holds either f _a,a or f _b,b read from the memory 81a, and outputs either f _a,a or f _b,b at the timing when the other of f _a,a or f _b,b is read from the memory 81a, and supplies it to the difference calculation circuit 83a.

The register 82b holds either _{dφ(a), φ(a)} or dφ _{(b), φ(b)} read from the memory 81b, and outputs either _{dφ(a), φ(a) or dφ(b) , φ(b)} at the timing when the other of _{dφ(a), φ(a)} or _{dφ(b), φ(b)} is read from the memory 81b, and supplies it to the difference calculation circuit 83b.

The difference calculation circuit 83a calculates f _b,b −f _a,a in equation (44).
The difference calculation circuit 83b calculates d _{φ(a), φ(a)} −d _{φ(b), φ(b)} in equation (44).

The multiplication circuit 84 calculates (f _b,b − _{f a,a} ) (d _{φ(a),φ(a) − d φ(b), φ(b) ) in equation (44) based on the calculation results of the difference calculation circuits 83 a and 83 b} .
Compensation term calculation circuit 85 calculates the compensation terms (terms that compensate for the elimination of the condition that the calculation of i = a, b is skipped), which are the four items on the right side of equation (44). To calculate the compensation terms, compensation term calculation circuit 85 receives f _a,b , d _{φ(a), and φ(b)} from multiplexers 33 a and 33 b in replica processing circuit 70 a 1, and receives f b _,a , d φ(b _{), and φ(a)} from multiplexers 33 a and 33 b in replica processing circuit 70 a 2.

Adder circuit 86 receives the value of one item on the right-hand side of equation (44) from dot product calculation circuit 34 of replica processing circuit 70a1, and receives the values of two items on the right-hand side of equation (44) from dot product calculation circuit 34 of replica processing circuit 70a2. Adder circuit 86 also receives the values of three items on the right-hand side of equation (44) from multiplier circuit 84, and the value of four items on the right-hand side of equation (44) from compensation term calculation circuit 85. Adder circuit 86 then calculates the sum of these values to generate and output ΔC _asym .

These circuits, memories, etc. may be included in the storage unit 11 or the processing unit 12 in FIG.
Note that the configuration for storing the integer allocation vector φ representing the allocation state, the configuration for determining whether or not to accept a proposal for an allocation change that will result in the calculated ΔC _ex , and the like are omitted from the illustration.

With the above hardware configuration, it is possible to perform a local search using the SAM method on a QAP. A configuration for performing a local search using the SAM method on a QSAP can also be realized by appropriately modifying the above hardware configuration. For example, an arithmetic circuit (such as a multiplication circuit or addition circuit) is added to perform the calculation on the second line of equation (24).

The above-described processing contents (for example, FIGS. 6 to 8, 10, 11, etc.) can be realized by software by causing the data processing device 10 to execute a program.
The program can be recorded on a computer-readable recording medium. Examples of recording media that can be used include magnetic disks, optical disks, magneto-optical disks, and semiconductor memories. Magnetic disks include flexible disks (FDs) and HDDs. Optical disks include compact discs (CDs), recordable/rewritable CD-Rs/RWs, digital versatile discs (DVDs), and DVD-Rs/RWs. The program may be recorded on a portable recording medium and distributed. In this case, the program may be copied from the portable recording medium to another recording medium and executed.

FIG. 31 is a diagram illustrating an example of hardware of a computer, which is an example of a data processing device.
The computer 90 includes a CPU 91, a RAM 92, a HDD 93, a GPU 94, an input interface 95, a media reader 96, and a communication interface 97. The above units are connected to a bus.

The CPU 91 is a processor including an arithmetic circuit that executes program instructions. The CPU 91 loads at least a portion of the program and data stored in the HDD 93 into the RAM 92 and executes the program. The CPU 91 may include multiple processor cores to execute the processing of multiple replicas in parallel, for example, as shown in FIG. 4. The computer 90 may also include multiple processors. A collection of multiple processors (a multiprocessor) may also be called a "processor."

RAM 92 is a volatile semiconductor memory that temporarily stores programs executed by CPU 91 and data used by CPU 91 for calculations. Note that computer 90 may be equipped with other types of memory besides RAM 92, or may be equipped with multiple memories.

The HDD 93 is a non-volatile storage device that stores software programs such as the OS (Operating System), middleware, and application software, as well as data. The programs include, for example, a program that causes the computer 90 to execute a process for searching for a solution to the assignment problem described above. The computer 90 may also be equipped with other types of storage devices, such as flash memory or an SSD (Solid State Drive), or may be equipped with multiple non-volatile storage devices.

In accordance with instructions from the CPU 91, the GPU 94 outputs an image (for example, an image showing the calculation results of an assignment problem) to a display 94a connected to the computer 90. The display 94a can be a CRT (Cathode Ray Tube) display, a liquid crystal display (LCD), a plasma display panel (PDP), an organic electro-luminescence (OEL) display, or the like.

The input interface 95 receives input signals from an input device 95a connected to the computer 90 and outputs them to the CPU 91. Examples of the input device 95a include pointing devices such as a mouse, touch panel, touchpad, or trackball, a keyboard, a remote controller, and button switches. Multiple types of input devices may also be connected to the computer 90.

The media reader 96 is a reading device that reads programs and data recorded on a recording medium 96a. Examples of recording media 96a that can be used include magnetic disks, optical disks, magneto-optical disks (MO: Magneto-Optical disks), and semiconductor memories. Magnetic disks include floppy disks and HDDs. Optical disks include CDs and DVDs.

The media reader 96 copies programs and data read from the recording medium 96a to other recording media such as the RAM 92 and HDD 93. The read programs are executed by the CPU 91, for example. The recording medium 96a may be a portable recording medium and may be used to distribute programs and data. The recording medium 96a and HDD 93 may also be referred to as computer-readable recording media.

The communication interface 97 is connected to the network 97a and communicates with other information processing devices via the network 97a. The communication interface 97 may be a wired communication interface connected to a communication device such as a switch via a cable, or a wireless communication interface connected to a base station via a wireless link.

The above describes one aspect of the program, data processing device, and data processing method of the present invention based on an embodiment, but these are merely examples and are not limited to the above description.

For example, in the above explanation, the columns of the distance matrix are rearranged according to the allocation status, but the same effect can be achieved by rearranging the rows of the distance matrix according to the allocation status, by modifying the formula appropriately.

10 Data processing device 11 Storage unit 12 Processing unit

Claims (9)

A program that causes a computer to execute a process of searching for a solution to an assignment problem by local search using an evaluation function that represents a cost according to an assignment state,
calculating, using vector arithmetic operations, a first change in the evaluation function when a first allocation change occurs in which allocation destinations of a first element and a second element among the plurality of elements are swapped, based on a flow matrix representing flow amounts between a plurality of elements allocated to a plurality of allocation destinations and a distance matrix representing distances between the plurality of allocation destinations, both of which are stored in a memory;
determining whether to allow the first allocation change based on the first change amount;
If it is determined that the first allocation change is permitted, the allocation state is updated, and the distance matrix is updated so that two columns or two rows corresponding to the first element and the second element are swapped.
A program that causes the computer to execute a process. The program of claim 1 causes the computer to execute a process that, when the assignment problem is QAP, determines whether to accept the first assignment change based on a comparison result between an acceptance probability calculated based on the first change amount and the value of a temperature parameter and a random value. Before calculating the first change amount,
When the assignment problem is QSAP, a second change in the evaluation function is calculated when a second assignment change occurs in which the first element is assigned to a first assignment destination;
determining whether to accept the second allocation change based on a comparison result between an acceptance probability calculated based on the second change amount and the value of the temperature parameter and a random value;
If it is determined that the second allocation change is permitted, the allocation state and the distance matrix are updated;
After calculating the first change amount,
determining whether or not to permit the first allocation change based on a comparison result between the first change amount and a predetermined value;
The program according to claim 1 , which causes the computer to execute a process. reading the distance matrix row by row from the memory;
Select two values in the two columns included in the read row;
writing the two values to the memory by swapping their storage locations;
Repeat the process to swap the two columns.
The program according to claim 1 , which causes the computer to execute a process. the memory further stores a transpose of the distance matrix;
storing a first row of the transposed matrix in a first shift register and a second row of the transposed matrix in a second shift register, the first row corresponding to the two columns of the distance matrix;
the storage locations of two values output from the first shift register and the second shift register are swapped and written into the memory, and the process is repeated to swap the two columns;
The program according to claim 1 , which causes the computer to execute a process. the memory includes a first memory and a second memory in which the distance matrix is stored;
reading the distance matrix row by row from the first memory;
Select two values in the two columns included in the read row;
writing the two values to the second memory by swapping their storage locations;
Repeat the process to swap the two columns.
The program according to claim 1 , which causes the computer to execute a process. the local search is performed using parallel tempering with a plurality of replicas, each replica having a different temperature parameter value;
the memory stores the distance matrices of a first replica and a second replica among the plurality of replicas;
calculating the first change amount based on the distance matrix of the second replica while updating the distance matrix of the first replica;
The program according to claim 1 , which causes the computer to execute a process. A data processing device that searches for a solution to an assignment problem by local search using an evaluation function that represents a cost according to an assignment state,
a storage unit that stores a flow matrix that represents flow amounts between a plurality of elements that are allocated to a plurality of allocation destinations, and a distance matrix that represents distances between the plurality of allocation destinations;
a processing unit that calculates, based on the flow matrix and the distance matrix, a first change amount of the evaluation function using vector arithmetic operation when a first allocation change occurs in which allocation destinations of a first element and a second element among the plurality of elements are swapped, determines whether or not to allow the first allocation change based on the first change amount, and when it is determined that the first allocation change is allowed, updates the allocation state and updates the distance matrix so that two columns or two rows corresponding to the first element and the second element are swapped;
A data processing device having: A data processing method for searching for a solution to an assignment problem by local search using an evaluation function that represents a cost according to an assignment state, comprising:
The computer
calculating, using vector arithmetic operations, a first change in the evaluation function when a first allocation change occurs in which allocation destinations of a first element and a second element among the plurality of elements are swapped, based on a flow matrix representing flow amounts between a plurality of elements allocated to a plurality of allocation destinations and a distance matrix representing distances between the plurality of allocation destinations, both of which are stored in a memory;
determining whether to allow the first allocation change based on the first change amount;
If it is determined that the first allocation change is permitted, the allocation state is updated, and the distance matrix is updated so that two columns or two rows corresponding to the first element and the second element are swapped.
Data processing methods.