JP2023001055A - Program, data processing apparatus, and data processing method - Google Patents

Abstract

To solve assignment problems at high speed.SOLUTION: A storage unit 11 stores a flow matrix representing flows between a plurality of entities to be assigned to a plurality of destinations, and a distance matrix representing distances between the plurality of destinations. A processing unit 12 calculates a first change in an evaluation function, which is to be caused by a first assignment change of exchanging the destinations of first and second entities among the plurality of entities, with vector arithmetic operations based on the flow and distance matrices, determines based on the first change whether to accept the first assignment change, and when determining to accept the first assignment change, updates an assignment state and updates the distance matrix by swapping the two columns or two rows (two columns in the example of FIG. 2) of the distance matrix corresponding to the first and second entities.SELECTED DRAWING: Figure 2

Description

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

Allocation problems are a class of NP-hard combinatorial optimization problems that have a variety of real-world applications such as vehicle dispatch routing and block placement in Field-Programmable Gate Arrays (FPGAs). An example of the assignment problem is the quadratic assignment problem (QAP) (see Non-Patent Document 1, for example).

In the secondary allocation problem, when n elements (facilities, etc.) are allocated to n allocation destinations, the amount of flow between elements (cost such as the amount of goods transported between facilities) and each element are allocated. This is a problem of finding an allocation that minimizes the sum of products of distances between allocation destinations. That is, the secondary assignment problem is a problem of searching for an assignment that minimizes the value of the evaluation function represented by Equation (1) below. The evaluation function expresses the cost according to the allocation state, and is also called a cost function.

In equation (1), f _i,j is the amount of flow between elements with identification numbers=i and j, and d _{φ(i) and φ(j)} are assigned to elements with identification numbers=i and j. Indicates the distance between assignees.
By the way, there is a Boltzmann machine (also called an Ising machine) using an Ising evaluation function as a device for calculating a large-scale discrete optimization problem, which is not good for von Neumann computers. A Boltzmann machine is a type of recurrent neural network.

A Boltzmann machine transforms a combinatorial optimization problem into an Ising model that represents the spin behavior of a magnetic material. Then, the Boltzmann machine uses a Markov chain Monte Carlo method such as pseudo-annealing or parallel tempering (see, for example, Non-Patent Document 2) to search for the state of the Ising model that minimizes the value of the Ising-type evaluation function (for example, , Non-Patent Document 3). The value of the evaluation function corresponds to energy. 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 represented by a combination of multiple state variable values (also called neuron values). 0 or 1 can be used as the value of each state variable.

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

One item on the right side is the value (0 or 1) of two state variables selected from N state variables without omission or duplication, and the weight value (of the two state variables) for all combinations of all state variables of the Ising model. It represents the strength of the interaction between s _i is the state variable with the identification number i, s _j is the state variable with the identification number j, and w _i,j indicates the magnitude of the interaction between the state variables with the identification numbers i and j. is a weight value. The two items on the right side are sums of products of bias coefficients and state variables for each identification number. b _i indicates the bias coefficient for identification number=i.

Also, the amount of change in energy (ΔE _i ) that accompanies the change in the value of s _i is represented by the following equation (3).

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

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

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

In equation (4), x is a vectorized state variable, 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 ). x _i,j =1 indicates that the element with identification number=i is assigned to the assignee with identification number=j, and x _i,j =0 indicates that the element with identification number=i = indicates that it is not assigned to the assignee of j.

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

U.S. 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

The method of calculating the quadratic assignment problem using the above evaluation function requires many iterations of processing such as calculation of the amount of change in energy and updating of the local field, and memory access is also repeated many times. time consuming.

An object of the present invention in one aspect is to provide a program, a data processing apparatus, and a data processing method capable of calculating allocation problems at high speed.

In one embodiment, a program for causing a computer to execute a process of searching for a solution to an allocation problem by local search using an evaluation function representing a cost according to an allocation state, wherein a plurality of allocations stored in a memory Based on a flow matrix representing the amount of flow between the plurality of elements to be allocated first and a distance matrix representing the distance between the plurality of allocation destinations, the first element and the second element among the plurality of elements calculating a first amount of change in the evaluation function using vector arithmetic operations when a first allocation change occurs in which an element allocation destination is replaced; and based on the first amount of change, the first If it is determined that the first allocation change is permitted, the allocation state is updated, and in the distance matrix, the first element and the second element A program is provided that causes the computer to perform a process of updating to swap two columns or two rows corresponding to .

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 enables fast computation of allocation problems.

It is a figure which shows the calculation example of QAP. FIG. 10 is a diagram showing an example of rearrangement of distance matrices and an example of a data processing device when calculating QAP; FIG. 10 is a diagram illustrating an example of cache matrix update calculation; FIG. 2 illustrates an example solver system that performs parallel tempering; FIG. 4 is a diagram showing an example of an algorithm for searching for a QAP solution by local search by parallel tempering; 4 is a flow chart showing the overall processing flow of local search by parallel tempering; FIG. 11 is a flow chart showing an example of the flow of replica initialization processing in the case of QAP; FIG. FIG. 11 is a flow chart showing an example of the flow of replica search processing in the case of QAP; FIG. FIG. 10 is a diagram showing an example of an algorithm for searching for a QSAP solution by local search by parallel tempering; FIG. 11 is a flow chart showing an example of the flow of replica initialization processing in the case of QSAP; FIG. FIG. 11 is a flow chart showing an example of the flow of replica search processing in the case of QSAP; FIG. FIG. 10 is a diagram showing evaluation results of the degree of speeding up of calculation by the SAM method and the BM$ method compared with the VΔC method; FIG. 10 is a diagram showing evaluation results of the degree of acceleration of vector-type arithmetic processing relative to scalar-type arithmetic processing; It is a figure which shows an example of load distribution. FIG. 10 is a diagram showing evaluation results of the degree of speeding up of arithmetic processing by load balancing; FIG. 4 is a diagram showing an example of a measurement algorithm; Fig. 2 shows the relative speedup measured for the VΔC, SAM, and BM$ methods and the memory hierarchy occupied as a function of problem size; FIG. 4 is a diagram showing an example of a ΔC generation circuit; FIG. 10 is a diagram illustrating a first example of a hardware configuration for permuting columns; It is a figure which shows the example of permuting a row|line|column. FIG. 10 is a diagram illustrating a second example of a hardware configuration for permuting columns; FIG. 10 is a diagram showing a first modification of the second example of the hardware configuration for permuting columns; FIG. 11 is a diagram showing a second modification of the second example of the hardware configuration for permuting columns; FIG. 13 is a diagram illustrating a third example of a hardware configuration for permuting columns; FIG. 12 is a diagram showing a modification of the third example of the hardware configuration for permuting columns; FIG. 10 is a diagram showing another example of a ΔC generation circuit; FIG. 10 is a diagram illustrating an example of a hardware configuration for processing two replicas; FIG. FIG. 10 is a diagram showing another example of a ΔC generation circuit; It is a figure which shows an example of a replica processing circuit. FIG. 4 is a diagram showing an example of a replica processing circuit used in QAP calculation using an asymmetric matrix; It is a figure which shows the hardware example of the computer which is an example of a data processing apparatus.

BEST MODE FOR CARRYING OUT THE INVENTION Hereinafter, embodiments for carrying out the invention will be described with reference to the drawings.
The data processing apparatus of the present embodiment searches for a solution of a quadratic assignment problem (QAP) or a quadratic semi-assignment problem (QSAP: Quadratic Semi-Assignment Problem) as an example of an assignment problem by local search. QAP, QSAP and local search are described below.

(QAP)
FIG. 1 is a diagram showing an example of QAP calculation. QAP, when assigning n elements (facilities, etc.) to n assignees, minimizes the sum of the products of the flow amount between the elements and the distance between the assignees to which each element is assigned. It is a matter of asking for allocation.

FIG. 1 shows an example of allocating four facilities with identification numbers 1 to 4 to four allocation destinations (L ₁ to L ₄ ).
A flow matrix representing the amount of flow between n elements is represented by the following equation (6).

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

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

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

QAP computation is done by searching for allocations that minimize equation (1) above. The allocation state of n elements to n allocation destinations is represented by an integer allocation vector φ or a binary state matrix X. φ is an element of the set Φ _n . Φ _n is the set of all permutations of the set N={1,2,3,...,n}. A binary variable x _i,j included in the binary state matrix X is represented by the following equation (8).

In FIG. 1, the facility with identification number = 1 is the allocation destination with identification number = 2, the facility with identification number = 2 is the allocation destination with identification number = 3, the facility with identification number = 3 is the allocation destination with identification number = 4, An example is shown in which facilities with identification number=4 are assigned to allocation destinations 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.

In QAP, one or both of the flow matrix and the distance matrix are symmetrical, and both the flow matrix and the distance matrix are asymmetrical. In this embodiment, the focus is mainly on QAP using a symmetric matrix (diagonal elements are 0 (biasless)). This is because such QAPs are the bulk of the instances and to simplify the calculations. However, QAP with symmetric matrices can be directly converted to QAP with asymmetric matrices.

(QSAP)
QSAP is a variant of QAP. In QSAP, the number of elements and the number of assignees are not equal. For example, it is permissible to assign multiple elements to each assignee. In QSAP, the distance matrix is represented by Equation (9) below.

The diagonal elements of the distance matrix (D) are non-zero to allow for intra-assignment routing. QSAP also uses an additional matrix B expressed in Equation (10) below.

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

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

(local search)
Local search searches for candidate solutions within neighboring states that are reachable by changing the current state. Pairwise exchange between elements is a method of performing local search for QAP. In this approach, two elements are selected and their assignments are exchanged. The value of the evaluation function (formula (1)) obtained by exchanging the allocation destinations (identification numbers = φ(a), φ(b)) of the two elements (identification numbers = a, b) can be expressed by the following equation (13).

As shown in equation (13), the amount of change (ΔC _ex ) is generated by a sum-of-products operation loop.
When local search is performed on 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 constant change in cost for allocating elements to each assignee, and is expressed by Equation (15) below.

The state space in QSAP is not constrained only to the aforementioned permutations representing allocation states. An element can be relocated to its assignee regardless of whether another element is assigned to that assignee. The amount of change in the value of the evaluation function when the element of identification number = a is changed from the current allocation destination (identification number = ψ(a)) to the allocation destination of identification number = l is given by the following formula (16 ), (17).

In the local search, based on the amount of change calculated as described above, it is determined whether or not to accept a proposal of assignment change that causes a change in the value of the evaluation function of the amount of change. Whether or not to accept the offer is determined based on predefined criteria, eg, greedy. If the greedy method is used, allocation changes that reduce the cost (value of the cost function (energy)) are accepted. When the greedy method is used, the Proposal Acceptance Rate (PAR) is high at the beginning of the search, but the search gets stuck at a local minimum of the evaluation function and no further improvement moves are found. It tends to approach 0 later.

Instead of the greedy method, this embodiment can use a probabilistic local search method such as the pseudo-annealing method. In the probabilistic local search method, the temperature parameter (T) is used to add randomness to the assignment changes, using the Metropolis method acceptance probability (P _acc ) expressed in equation (18) below. be able to.

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 at a local minimum of the evaluation function.

The calculation of ΔC for QAP and QSAP accounts for most of the total processing time in the data processing system as it searches for a solution to the allocation problem.
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. Calculation of ΔC _ex of QAP by Equation (13) can be easily vectorized using inner product calculation, provided that there is no condition to skip the calculation of i=a,b.

The data processing apparatus of the present embodiment eliminates this condition, calculates ΔC _ex as a product-sum operation loop for all n elements, and uses the flow amount and the distance as shown in the following equation (19). Compensate by removing that condition by adding an additional product of as a compensation term.

Equation (19) expresses the flow amount difference vector ΔF (equation (20)) and the distance difference vector ΔD (equation (21)) when the allocation destinations of the elements with identification numbers = a and b are switched. can be reformulated using inner products as in equation (22).

Δ ^b _a F in equation (20) represents the difference vector between rows b and a of the flow matrix, F _b,* denotes row b of the flow matrix, and F _a,* denotes 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 allocation state. This can be done mathematically by multiplying the transposed binary state matrix X (denoted X ^T ) by (D _φ(a),* −D _φ(b),* ) as in equation (21). shown D _φ(a),* denotes row a of the distance matrix and D _φ(b),* denotes row b of the distance matrix.

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

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

Since QSAP does not restrict multiple elements to be assigned to the same assignee, an additional compensation term is added to ΔC as in the second line of equation (24).
Also, the vectorized form of the calculation of ΔC ^s _rel by rearrangement corresponding to Equation (17) can be expressed as Equation (25) below.

Since the movement of the assignment destination of the element due to rearrangement is performed for one element, it is not necessary to calculate the vector ΔF and compensate by adding the product of the flow amount and the distance.
(Rearrangement of distance matrix)
The above VΔ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, the rearrangement of the elements of ΔD performed by the VΔC method can be executed in a time proportional to the number of elements = n, but it is executed each time the calculation of ΔC is performed by switching the assignment destination of the two elements. Poor processing efficiency. It can also be computationally expensive.

In order to minimize the computational cost of generating ΔD with proper element order, the data processing apparatus of the present embodiment rearranges the columns of the distance matrix according to the current allocation state. Hereinafter, this method may be called SAM (State-Aligned D Matrix) method. Also, hereinafter, the distance matrix rearranged in this way is called a state-ordered D matrix, and is denoted as D ^X when the allocation problem is QAP and D ^S when the allocation problem is QSAP.

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

The storage unit 11 is, for example, a volatile storage device that is an electronic circuit such as a DRAM (Dynamic Random Access Memory), or a non-volatile storage device that is an electronic circuit such as a HDD (Hard Disk Drive) or flash memory. . The storage unit 11 may 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 and the like, and also stores the aforementioned flow matrix, distance matrix, state-aligned D matrix (D ^X or D ^S ), allocation state (integer allocation vector φ or binary state matrix X).

The processing unit 12 can be realized by a processor, which is hardware such as a CPU, GPU, DSP (Digital Signal Processor), for example. Also, the processing unit 12 may be implemented 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. Note that the processing unit 12 may be a set of multiple processors.

The processing unit 12 repeats the process of interchanging the allocation destinations of the two elements, and searches for an allocation state in which the value of the evaluation function shown in Equation (1) or Equation (11) is minimized, for example. The optimum solution is the allocation state that has the minimum value among the minimum values of the evaluation function. 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 optimal solution).

FIG. 2 shows an example of permutation of the distance matrix when calculating the QAP. An allocation state for a certain time step=t and a state-aligned D matrix having an arrangement corresponding to the allocation state are denoted by φ ^(t) and D ^X(t) . If an allocation change is proposed and accepted to replace the allocation destinations of the elements with identification numbers=1 and 3, at time step=t+1, the allocation state is as represented by φ ^(t+1) in FIG. A state-aligned D-matrix is then generated, ^D_X(t+1) , with the first and third columns swapped with respect to ^D_X(t) to correspond to the accepted allocation states.

Note that if the proposed allocation change is not accepted, φ ^(t) and D ^X(t) are not updated.
For QAP, the state-aligned D matrix (D ^X ) is the transposed matrix X ^T of the binary state matrix X representing the current allocation state, and the original distance matrix (D), as shown in equation (26) below. can be expressed by multiplying An example of a hardware configuration in which the columns of the state-aligned D-matrix are permuted in terms of hardware will be described later.

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

For QSAP, the state-aligned D matrix (D ^S ) is the transposed matrix S ^T of the binary state matrix S representing the current allocation state, and the original distance matrix (D ), as shown in equation (28) below. can be expressed by multiplying D ^S is a matrix with m rows and n columns.

The vector ΔD ^S due to the allocation target exchange and rearrangement can be calculated by the following equations (29) and (30), respectively, without rearranging the elements of ΔD, and the above equations (24) and ( 25), ΔC ^s _ex and ΔC ^s _rel can be calculated.

In the VΔC method described above, to rearrange the elements of the vector ΔD, an operation proportional to the number of elements=n is performed each time ΔC is calculated. In contrast, in the above approach (the SAM method), D ^X and D ^S are reordered to match the allocation state only if the proposed allocation change is accepted. . This increases computational efficiency by reducing the number of operations required per iteration on average, thus reducing computation time and allowing the allocation problem to be computed quickly.

(Boltzmann machine caching method (comparative example))
The following method is also conceivable as a method of vectorizing the ΔC calculation. This technique is a method of calculating and storing partial ΔC values (corresponding to the local fields (h _i ) shown in equation (3) above) corresponding to bit inversions of the binary state matrix. Hereinafter, this method may be called the Boltzmann machine caching method (abbreviated as BM$ method).

A cached local field is used for each bit of the binary state matrix (X or S). When the allocation change proposal is accepted, the matrix with n rows and ⁿ columns by the local field (hereinafter referred to as the cache matrix) is updated in a time proportional to n2. can be done in a time independent of

For QAP ^, the cache matrix (H) is generated by equation (31) at the start of the search process in time proportional to n3.

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

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

FIG. 3 is a diagram illustrating an example of cache matrix update calculation.
Updates are made one row of the cache matrix at a time. Processing is skipped for rows corresponding to 0 elements of vector ΔF.

For QSAP, the cache matrix (H ^s ) is generated by equation (34) at the start of the search process.

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

If the reassignment proposal is accepted, the cache matrix is updated by 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 according to equation (37). Then, ΔC ^s _rel can be calculated by substituting the generated ΔF·ΔD into the above equation (25).

If the relocation assignment change proposal is accepted, the cache matrix is updated according to equation (38).

In the BM$ method as described above, a cache matrix is held, and as the scale of the problem increases, the capacity of memory for holding the cache matrix also increases. This makes it difficult to use a memory with a relatively small capacity that can be read out at high speed. The SAM method need not maintain 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, BM$ method), a solver system as shown below was designed. In the following description, it is assumed that the processing unit 12 of the data processing device 10 includes a multi-core CPU having SIMD functions, and a parallel tempering algorithm is implemented on the multi-core CPU.

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

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

In the case of QAP, for each replica, the aforementioned integer allocation vector φ, cache matrix (H) (for BM$ method), state-aligned D matrix (D ^X ) (for SAM method), evaluation function value (C ) is retained.

In the case of QSAP, for each replica, the above-mentioned integer assignment vector ψ, cache matrix (H ^s ) (for BM$ method), state-aligned D matrix (D ^s ) (for SAM method), evaluation function value ( C) is retained.

The flow matrix (F) shown in Equation (6) and the matrix B shown in Equation (10) used for QSAP are common to all replicas 20b1 to 20bM.
Different temperature parameters (T) are set for the M replicas 20b1 to 20bM. The replicas 20b1 to 20bM are set to any temperature that increases stepwise from T _min to T _max . Hereinafter, such a stepwise increasing temperature parameter value may be referred to as a temperature ladder.

The parallel tempering controller 21 performs the following operations based on the evaluation function values (C _k , C _k+1 ) and T _k , T _k+1 between the adjacent replicas for which the temperature parameter values (T _k , T _k+1 ) are set. The values of T _k and T _k+1 are exchanged with the exchange admissible probability SAP expressed by Equation (39).

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

FIG. 5 is a diagram showing an example of an algorithm for searching a QAP solution by local search by parallel tempering. Note that in the example of the algorithm in FIG. 5, the number of replicas for which local searches are performed in parallel is 32. In FIG. 5, "(1)", "(31)", etc. shown in the 16th line, the 18th line, etc. represent the aforementioned formulas (1), (31), and the like.

The algorithm begins by initializing the values of the temperature parameters and the replica state. The solver system then enters an optimization loop where I iterations of the local search are performed on all replicas in parallel.

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

In the QAP local search, the loop is run for a preset number of iterations (I), two elements (facility in the example of FIG. 5) are selected, and ΔC _ex is calculated by the chosen method (VΔC method, BM$ method or SAM method).

P _acc in equation (18) is then calculated using the calculated ΔC _ex and the value of the current temperature parameter for each replica. A random value in the range 0 to 1 is then generated and compared to the P _acc and a Bernoulli trial is performed to determine if the proposal for swapping assignments between the two selected elements is acceptable. . If the offer is accepted, the status (assignment status) is updated.

FIG. 6 is a flow chart showing the overall processing flow of local search by parallel tempering. FIG. 6 shows the overall processing flow of local search by parallel tempering of the algorithm shown in FIG.

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

Next, an optimization loop (steps S14 to S16) is performed while increasing i by one between i=0 and i<M−1.
In the optimization loop, replica search processing is performed (step S15). In the replica search process, local searches are performed I times for each replica using the set temperature parameter value.

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

Note that the data processing device 10 may output search results (for example, minimum values, C _min , φ _min described later) when the search ends. The search result may be displayed on a display device connected to the data processing device 10, for example, or may be transmitted to an external device.

Although FIG. 6 shows an example in which processing is performed sequentially for each replica, for example, a processor having 32 cores can execute the processing of steps S11, S12, and S15 for 32 replicas in parallel. be.

Next, the flow of the replica initialization process in step S12 and the replica search process in step S15 will be described using a flow chart, taking the case of executing the SAM method as an example.
FIG. 7 is a flowchart showing an example of the flow of replica initialization processing in the case of QAP.

First, the integer assignment vector φ is initialized (step S20). In the processing of step S20, for example, the allocation destination of the element (placement destination of the facility) is determined at random. Then, the initial cost (C) is calculated by the calculation of formula (1) (step S21).

After that, the initialized φ and C initialize the minimum values (C _min and φ _min ) (step S22). Also, the state-aligned D-matrix (D ^X ) is initialized so as to correspond to the initial value of φ (step S23). After that, the process returns to the process of 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 between i=1 and 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 exchanged is calculated by Equation (22) (step S32).

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

If it is determined that P _acc >rand( ), then the permutation of columns a and b of D ^X (D ^X _{*, a} , D ^X _{*, b} ←D ^X _{*, b} , D ^X _{*, a} ) is is performed (step S35). Furthermore, the allocation state is updated (φ(a), φ(b)←φ(b), φ(a)) by replacing the allocation destination (step S36), and the evaluation function value (C) is updated ( C←C+ΔC _ex (a, b)) is performed (step S37).

After that, it is determined whether or not 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 determined, or if it is determined that C<C _min is not determined in the process of step S38, until i=I The processing from step S31 is repeated. When i=I, the process returns to the flowchart shown in FIG.

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

The local search function that computes the QSAP shown in FIG. and before the exchange loop that does the update.

Next, the flow of the replica initialization process in step S12 and the replica search process in step S15 in FIG. 6 when calculating QSAP will be described with reference to a flowchart, taking the case of executing the SAM method as an example. .

FIG. 10 is a flowchart showing an example of the flow of replica initialization processing in the case of QSAP.
First, the integer assignment vector ψ is initialized (step S50). In the processing of step S50, for example, the allocation destination of the element (location of the facility) is determined at random. Then, the initial cost (C) is calculated according to the equation (11) (step S51).

After that, the minimum values (C _min and ψ _min ) are initialized by the initialized ψ and C (step S52). Also, the state-aligned D-matrix (D ^S ) is initialized according to equation (28) so as to correspond to the initial value of ψ (step S53). After that, the process returns to the process of 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.
Also in the replica search process in the case of QSAP, an iteration loop (steps S60 to S78) is performed while increasing i by one between i=1 and i<I.

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

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

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

After that, it is determined whether or not 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 when it is determined in the process of step S64 that P _acc >rand( ) is not obtained, or when it is determined in the process of step S68 that C<C _min is not satisfied, the process of step S70 is performed. will be

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

Next, it is determined whether or not ΔC ^s _ex <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, the permutation of columns a and b of D ^S (D ^S _{*, a} , D ^S _{*, b} ← D ^S _{*, b} , D ^S _{*, a} ) is a row. (step S73). Furthermore, the allocation state is updated (ψ(a), ψ(b)←ψ(b), ψ(a)) by replacing the allocation destination (step S74), and the evaluation function value (C) is updated ( C←C+ΔC ^s _ex ) is performed (step S75).

After that, it is determined whether or not C<C _min (step S76), and if it is determined that C<C _min , the minimum value is updated (C _min ←C, ψ _min ←ψ) ( step S77).

After the process of step S77, or if it is determined that ΔC ^s _ex <0 is not determined in the process of step S72, or if it is determined that C<C _min is not determined in the process of step S76, steps are performed until i=I. The processing from S61 is repeated. When i=I, the process returns to the flowchart shown in FIG.

The order of processing shown in FIGS. 10 and 11 is an example, and the order of processing may be changed as appropriate.
(Evaluation result of calculation speed)
First, the results of evaluating the calculation speed by three methods, the VΔC method, the SAM method, and the BM$ method, using scalar arithmetic processing, are shown. Calculation targets were 10 QAP instances, and the calculation was performed by the algorithm shown in FIG. The calculations were performed by parallel tempering of 64 replicas using two solver systems as shown in FIG.

Each instance was run sequentially 100 times with an independent random number seed and the time to reach BKS (TtS) was recorded.
FIG. 12 is a diagram showing evaluation results of the degree of speeding up of calculation by the SAM method and the BM$ method compared with 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 BM$ method compared to the VΔC method. there is The degree of speedup is represented by the ratio of TtS by the SAM method and BM$ method to TtS by the VΔC method.

With the SAM method, the computational speed is faster than with the ^VΔC method, despite the additional computational cost of updating DX.
When the BM$ method is used, the calculation speed is improved by a factor of 2.57 in terms of geometric mean as compared with the VΔC method.

(Improved calculation speed in case of SIMD)
Next, the result of evaluation of the degree of improvement in computational speed of vector-type operation processing relative to scalar-type operation processing will be shown. SIMD built-in functions of AVX (Advanced Vector eXtensions) 2 were used to perform vector type arithmetic processing.

The same 10 QAP instances were run with vector-based computation using a procedure similar to that described above.
FIG. 13 is a diagram showing an evaluation result of the degree of acceleration 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 acceleration of vector-type operation processing relative to scalar-type operation processing. there is The degree of speedup is represented by the ratio of TtS due to vector-type arithmetic processing to TtS due to scalar-type arithmetic processing for each of the VΔC method, SAM method, and BM$ method.

On average, the VΔC method was twice as fast as the scalar type arithmetic processing, and the BM$ method and the SAM method were three times or more faster than the scalar type arithmetic processing. The V.DELTA.C method has the least advantage of SIMD built-in functions because the reordering of the elements of .DELTA.D performed in the V.DELTA.C method occupies a significant portion of the execution time as described above.

(dynamic load balancing)
In parallel tempering, the PAR differs between the replica set at the lowest temperature and the replica set at the highest temperature. The higher the value of T (higher the temperature), the higher the PAR. This can cause the SAM and BM$ update processes to cause large execution time gaps between the threads that process the replicas.

To mitigate this, data processor 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.

FIG. 14 is a diagram illustrating an example of load distribution. FIG. 14 shows calculation time of ΔC before exchange of temperature parameters and update processing time for 32 replicas for which temperature ladders of T1 to T32 (T1<T2< . . . <T32) are set.

When load balancing is not performed, replicas with larger temperature parameter values take longer to update, and replicas with smaller temperature parameters take less time to update. For this reason, replicas with smaller temperature parameter values have longer idle times.

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

As a result, as shown in FIG. 14, the execution times of the threads can be made substantially equal among the replicas, and the idle time can be shortened.
In order to quantitatively demonstrate the effectiveness of such load distribution, the above-mentioned 10 instances were used as calculation targets, and the results of evaluating the degree of speed-up of arithmetic processing due to load distribution are shown below.

FIG. 15 is a diagram showing evaluation results of the degree of speeding up of arithmetic processing by load distribution. The horizontal axis represents 10 instances of QAP and the geometric mean "GEOMEAN", and the vertical axis represents the degree of speed-up when load balancing is performed as compared to the case when load balancing is not performed. The degree of speedup is determined by the ratio of TtS when the SAM method and BM$ method are executed with load distribution to the TtS when the SAM method and BM$ method are executed without load distribution. is represented. Note that the same set of temperature parameter values (for example, T1 to T32 in FIG. 14) are used both with and without load sharing.

It was found that by load distribution, an average speed-up of 1.86 times was achieved in the SAM method and an average speed-up of 1.38 times was achieved in the BM$ method. Differences in load balancing effectiveness between instance families may be due to different PAR gaps between extremes of temperature parameters due to instance family characteristics.

Table 1 summarizes the results of evaluating the computation speedup provided by each of the above functions.

Two degrees of speedup (Incremental Speed-Up and Cumulative Speed-Up) are shown. The degree of speedup was calculated using the geometric mean of TtS for the 10 instances described above. there is

(Benchmark result)
With respect to the VΔC method, the SAM method, and the BM$ method, benchmark scores were measured using a data processing apparatus 10 having a predetermined hardware configuration including a 64-core CPU.

An instance in the QAP library reference (see Non-Patent Document 5) was used to measure the benchmark score of the QAP solver performing the above three methods. Furthermore, the instances proposed in Non-Patent Documents 6 and 7 were used. Also, in the present embodiment, an attempt was made to find a new BKS state for a large problem size QAP (no known optimal solution) given in Non-Patent Document 7. The set of instances introduced in Non-Patent Document 8 was also used to measure the benchmark score of the QSAP solver.

(QAP benchmark)

Table 2 shows the TtS for each instance of QAP by solvers implementing three methods of vector-based computation (VΔC method, SAM method with load balancing, and BM$ method). Furthermore, Table 2 shows the degree of speed-up (calculated using the geometric mean of TtS) with TtS as the reference (1.00) by vector type arithmetic processing using the VΔC method. Furthermore, Table 2 shows the TtS and the degree of speedup by two open solvers, ParEOTS (see Non-Patent Document 9) and PBM (Permutational Boltzmann Machine) (see Non-Patent Document 1). there is These two published solvers are capable of solving some difficult instances with 100% success rate within a 5 minute timeout window.

The TtS values shown in Table 2 are the mean TtS values (with a 99% confidence interval of value). The timeout limit for each run is 5 minutes. The TtS values for ParEOTS and PBM are according to their respective disclosures.

As shown in Table 2, the BM$ method showed the best performance across all instances among the five solvers, speeding up by 2x over PBM and over 300x over ParEOTS. showed that. The SAM and BM$ methods showed an average speed improvement of 1.92 and 3.22 times over the VΔC method, respectively. Note that the difference in results in Table 2 compared to Table 1 is due to the instance used.

However, which of the SAM method and the 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 described later (see FIG. 17), the SAM method tends to be superior to the BM$ method when the PAR is high. This evaluation result is the result of mounting on a CPU, and is the result of a system having a relatively large memory. In the case of implementation in a dedicated circuit, the SAM method, which does not require storing local fields, may be superior to the BM$ method.

(QSAP benchmark)
Compared to QAP, there are few previous disclosures on QSAP solvers. The main reference material is the PMITS (Parallel Memetic Iterative Tabu Search) algorithm implemented in a CPU with 20 cores (see Non-Patent Document 8). PMITS uses 50% of the population to arrive at the same best solution as the stopping criterion with a timeout limit of 1 hour. This method is a reasonable way to predict convergence in the Memetic algorithm using cooperative replicas. However, in parallel tempering, where replicas are in constant flow along a temperature ladder, such criteria are not valid.

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

As shown in Table 3, when compared to QAP, the performance of the SAM and BM$ methods compared to the VΔC method across QSAP instances is nearly doubled. This is believed to be due to the higher PAR.

(QAP scaling on extended Taillard instances)
Some QAP instances, introduced in Non-Patent Document 7, have unknown optimal solutions and are of size up to n=729. Rarely used for benchmarking because of its size and difficulty, data processor 10 uses the VΔC, SAM, and BM$ methods within a given time limit to find a better solution. and ran these instances.

Previous attempts to improve the BKS values for these instances yielded little or no improvement despite running for minutes to hours (see, eg, Non-Patent Document 10).

In the experiments, each instance was run 20 times and finished 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 together with the best cost (value of the evaluation function of equation (1)) when applying the VΔC method, the SAM method and the BM$ method for each instance.

It was found that despite the short run times, the BKS for all but four instances of the SAM and BM$ methods could be improved, and the average cost was also lower than the previous BKS.

(scaling analysis)
Based on benchmark results and a qualitative comparison between the VΔC method and the SAM method, the SAM method is believed to be more efficient. The VΔC method continuously rearranges the elements of ΔD for each iteration, whereas the SAM method rearranges D ^X and D ^S to match the allocation state. This is because it is the only case where the assigned assignment change is accepted.

On the other hand, the relative performance of SAM and BM$ methods mainly depends on PAR. PAR can vary greatly depending on the search algorithm used. Moreover, even within the same search algorithm, the PAR can change at run-time using algorithms such as simulated annealing and between concurrently searched instances such as parallel tempering.

The results of comparison of execution times according to problem sizes and PAR values are shown below.
FIG. 16 is a diagram showing an example of a measurement algorithm.
The measurement algorithm first generates a flow matrix (F) and a distance matrix (D) with random values from 1-100. Then the allocation state (φ) is initialized. Next, processing is executed for a predetermined number of iterations (I _limit ), and the execution time of one loop is measured. Proposals are then randomly accepted at the desired PAR value.

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 of 100M and 10M are used, respectively. For problem sizes of n>1024, an _{I_limit} of 1M is used.

At each iteration, a value of ΔC _ex is calculated and used by a random number generator to determine whether or not to accept the proposal with the desired PAR.
64 loop instances were executed in parallel (one per core) to measure the performance while the CPU pipeline, cache, and memory were under full load. 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) of the BM$ method, two separate sets of simulations were performed. One is a fully dense flow matrix and the other is a sparse flow matrix with only 10% non-zero values. Each measurement was calculated for a total of 290 parameter combinations, combining 10 problem sizes in the range n=[100,5000] and 29 PARs in the range PAR=[0.0001,0.1]. was done.

FIG. 17 shows the relative speedups measured for the VΔC method, the SAM method, and the BM$ method, and the memory hierarchy occupied according to the problem size. In the diagram showing the degree of speedup, the horizontal axis represents PAR [%] and the vertical axis represents the degree of speedup (Speed-Up). FIG. 17 shows the degree of acceleration of the BM$ method with respect to the SAM method and the VΔC method when flow matrices with different densities are used. Furthermore, FIG. 17 shows the degree of speedup of the SAM method relative to the VΔC method. In the example of FIG. 17, the memory hierarchy includes L2 cache, L3 cache, and DRAM in descending order of storage capacity.

It should be noted that the maximum PAR values for a given flow matrix density from non-load-sharing parallel tempering simulations for 10 instances of QAP were as shown in Table 6 below.

Note that the QSAP simulation is omitted because it is similar to the QAP simulation. It is considered that there is no significant difference between the results of both simulations.
The relative speedup with problem size can be divided into three groups based on which tier of the memory hierarchy contains most of the data used for the search. Each core has a dedicated L2 cache (256 kB storage capacity in this calculation example), and a group of four cores share an L3 cache (16 MB storage capacity in this calculation example).

As shown in FIG. 17, in the case of the SAM method and the BM$ method, the ^DX and cache matrix (H) of each loop thread fit in the L2 cache of the core with a maximum size of n=256, and the problem up to n=1024 is fits in the L3 cache.

As shown in Figure 17, the relative performance (degree of speedup) between the BM$ method and the other two methods decreases as the problem size increases and memory hierarchies with larger storage capacities are used. do. When the data used for searching is moved to a higher memory hierarchy, the value of PAR required for the BM$ method to maintain superior performance over the SAM and VΔC methods is greatly reduced.

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

The speedup of the SAM method over the VΔC method can be divided into two groups based on problem size. If the problem size is ^n≤800 and most of the DX per loop instance fits in the cache, the SAM method has a clear advantage over VΔC in the range of PAR<10%. Moving the data used for lookups from the L2 cache to the L3 cache does not significantly affect the relative performance between the two methods.

(Hardware example)
An example of hardware for implementing the SAM method will be described below.
To simplify the description, the following description assumes that both the flow matrix and the distance matrix are symmetric matrices (diagonal elements are 0 (biasless)). As mentioned above, such QAPs are the bulk of the instances and to simplify the calculations. QAP with symmetric matrices can be directly converted to QAP with asymmetric matrices.

A QAP evaluation function using only a symmetric matrix can be expressed as in Equation (40) below.

The only difference between formula (1) and formula (40) is that "j=1" in formula (1) is "j=i".
When the evaluation function represented by Equation (40) is used, ΔC _ex used for QAP calculation can be represented by the following Equations (41) and (42) instead of Equations (19) and (22). .

FIG. 18 is a diagram showing an example of a ΔC generation circuit.
The ΔC generation circuit 30 has a flow matrix memory 31 a , a state-aligned 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 and memories included in the storage unit 11 or the processing unit 12 shown in FIG.

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

In the example of FIG. 18, the flow matrix memory 31a and the state alignment D matrix memory 31b are dual port memories having two ports. Such dual port memories are available 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 inner product calculation circuit 34 calculates the inner product of ^{Δ b} _a F and Δ ^φ(a) _{φ(b) DX} . The inner product calculation circuit 34 can be implemented, for example, by a plurality of multipliers connected in parallel.

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

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, a permutation (assignment change) of the assignment destination of the element that generates ΔC _ex is proposed, and if accepted, the state-aligned D matrix is such that the columns correspond to the accepted assignment state. be replaced.

FIG. 19 is a diagram showing a first example of a hardware configuration for permuting columns.
As shown in FIG. 19, the permutation of the columns of the state-aligned D-matrix stored in the state-aligned D-matrix memory 31b can be performed using multiplexers 40a and 40b and a switch 41, for example. These circuits may also be included in the processing unit 12 of FIG.

The multiplexer 40a sequentially selects and outputs the values of the first columns of the rows of the state-aligned D-matrix read out for permutation from the state-aligned D-matrix memory 31b.
The multiplexer 40b sequentially selects and outputs the values of the second columns of the rows of the state-aligned D-matrix read from the state-aligned D-matrix memory 31b.

The switch 41 writes the value output from the multiplexer 40a to the location where the value output from the multiplexer 40b was stored in the state-arranged D-matrix memory 31b. Also, the switch 41 writes the value output from the multiplexer 40b to the location where the value output from the multiplexer 40a was stored in the state-arranged D-matrix memory 31b.

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

First, d _1,3 and d _1,1 are selected by multiplexers 40 a and 40 b , and switch 41 switches the memory location. Next, d _2,3 and d _2,1 are selected by multiplexers 40 a and 40 b and switched by switch 41 . The same process is repeated for a total of n cycles to complete the column permutation.

FIG. 21 is a diagram showing a second example of a hardware configuration for permuting columns. In FIG. 21, the same reference numerals are assigned to the same elements as those shown in FIG.
As shown in FIG. 21, the state-aligned D-matrix memory 31b stores not only the state-aligned D-matrix but also the transposed matrix ((D ^X ) ^T ) of the state-aligned D-matrix. The permutation of the columns of the state-arranged D-matrix can be performed by using the elements of the read-out transposed matrix, for example, by the shift registers 45a and 45b in addition to the switch 41 described above. 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 in which the state-aligned D-matrix is permuted are read from the state-aligned 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-aligned D-matrix memory 31b, shifts them one cycle at a time, and outputs the values one at a time.

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

The switch 41 switches the storage location in the state-aligned D-matrix memory 31b so that the value output from the shift register 45a is stored in the location where the value output from the shift register 45b was stored. Also, the switch 41 switches the storage location in the state-arranged D-matrix memory 31b so that the value output from the shift register 45b is stored in the location where the value output from the shift register 45a was stored.

Even with such a hardware configuration, columns can be exchanged in n cycles. Also, in the case of the hardware configuration of FIG. 21, there is a possibility that wiring performance is improved compared to the hardware configuration of FIG. Moreover, the multiplexers 40a and 40b as shown in FIG. 19, which may have a relatively large circuit scale, are not required.

FIG. 22 is a diagram showing a first modification of the second example of the hardware configuration for permuting columns. In FIG. 22, the same symbols are assigned to the same elements as those shown in FIG.

In the hardware configuration of FIG. 22, the transposed matrix ((D ^X ) ^T ) of the state-aligned D matrix is stored in the transposed matrix memory 46 separated from the state-aligned D matrix memory 31b. configuration is different. From the transposed matrix memory 46, the above two rows of the transposed matrix are read. Other configurations are the same as in FIG.

FIG. 23 is a diagram showing a second modification of the second example of the hardware configuration for permuting columns. In FIG. 23, the same reference numerals are assigned to the same elements as those shown in FIG.

In the hardware configuration of FIG. 23, the flow matrix memory 31a shown in FIG. 18 stores not only the flow matrix but also the transposed matrix ((D ^X ) ^T ) of the state-aligned D matrix. The permutation of the columns of the state-aligned D-matrix can be performed by the switch 41 and the 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 permuting columns. In FIG. 24, the same reference numerals are assigned to the same elements as those shown in FIG.
In the hardware configuration of FIG. 24, the state-aligned D-matrices are stored in two state-aligned D-matrix memories 31b1 and 31b2.

Using the elements of each row of the state-aligned D matrix read from the state-aligned D matrix memory 31b1, the hardware configuration of FIG. 19 is used for the state-aligned D matrix stored in the state-aligned D matrix memory 31b2. The columns are exchanged as if

The state-aligned D-matrix after the permutation of columns is copied to the state-aligned D-matrix memory 31b1.
Even with such a hardware configuration, columns can be exchanged in n cycles. Also, in the case of the hardware configuration of FIG. 24, there is a possibility that wiring performance is improved compared to the hardware configuration of FIG.

FIG. 25 is a diagram showing a modification of the third example of the hardware configuration for permuting columns. In FIG. 25, the same reference numerals are given to the same elements as those shown in FIG.
In the hardware configuration of FIG. 25, the flow matrix memory 31a shown in FIG. 18 stores not only the flow matrix but also the transposed matrix ((D ^X ) ^T ) of the state-aligned D matrix. The permutation of the columns of the state-aligned D-matrix can be performed by the multiplexers 40a, 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 the ΔC generation circuit. In FIG. 26, the same reference numerals are given to the same elements as those shown in FIG.
In the ΔC generation circuit 50 of FIG. 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 which holds Fa _,* , outputs Fa,* at the timing when Fb _,* is read from the flow matrix memory 51a, and supplies it to the difference calculation circuit 32a and the multiplexer 33a. Also, D ^X _φ(a),* is held, D ^X _φ(a),* is output at the timing when D ^X _φ(b),* is read from the state-aligned D matrix memory 51b, and the difference calculation circuit 32b and a register 52b that feeds the multiplexer 33b.

Other configurations are the same as those in FIG.
By the way, 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. becomes possible.

FIG. 27 is a diagram showing an example of a hardware configuration for processing two replicas. FIG. 27 shows a configuration in which the .DELTA.C generation circuit 50 shown in FIG. 26 is combined with the hardware configuration for column exchange shown in FIG.

In the example of FIG. 27, the state-aligned D matrix memory 51b stores state-aligned D matrices (D ^X _R1 , D ^X _R2 ) of two replicas (R1, R2).
Allocation changes are accepted at one replica, and column swaps are performed for one replica using the write port of the state-aligned D-matrix memory 51b. Meanwhile, for the other replica, processing is performed to calculate ΔC _ex based on the elements of the state-aligned D-matrix read from the read port of the state-aligned D-matrix memory 51b.

Note that if neither replica is executing an update process, then the process for one replica is stalled while the ΔC _ex calculation is being performed for the other replica.

FIG. 28 is a diagram showing another example of the ΔC generation circuit. In FIG. 28, the same reference numerals are given to the same elements as those shown in FIG.
If only the symmetrical permutation problem is to be calculated and ΔC _ex is represented by the following equation (43), a ΔC generating circuit 60 as shown in FIG. 28 can also be used.

The ΔC generation circuit 60 does not require the multiplexers 33a and 33b as shown in FIG. 18, and has a selector 61 instead.
For example, the difference calculation circuit 32b has n subtractors 32bi for calculating d _{φ(a), i} −d _{φ(b), i} , and the selector 61 selects one of the output of the subtractor 32bi and 0. It has n 2-input 1-output multiplexers 61ai that select and output.

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

(Replica processing circuit)
A circuit that performs processing for each replica is a combination of any of the ΔC generation circuits 40, 50, and 60 described above and any of the hardware configurations described above for permuting the columns of the state-aligned D matrix. realizable.

FIG. 29 is a diagram illustrating an example of a replica processing circuit; FIG. 29 shows an example of a replica processing circuit 70 in which the hardware configuration for column exchange shown in FIG. 19 is combined with the ΔC generation circuit 50 shown in FIG. In FIG. 29, the same elements as those shown in FIGS. 19 and 26 are denoted by the same reference numerals.

In the example of FIG. 29, multiplexer 33b also has the function of multiplexer 40a shown in FIG.
These circuits can be included in the storage unit 11 or processing unit 12 of 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 the allocation change proposal that causes the calculated ΔC _ex are omitted from the drawing.

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

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

The replica processing circuit 80 for calculating ΔC _asym has two replica processing circuits 70a1 and 70a2 for calculating ΔC _ex shown in FIG. However, the flow matrix memory 51a of one replica processing circuit 70a2 stores the transposed matrix (F ^T ) of the flow matrix, and the state-aligned D matrix memory 51b of the replica processing circuit 70a2 stores the state-aligned D matrix A transposed matrix ((D ^X ) ^T ) is stored.

The replica processing circuit 80 further has 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.

Memory 81a stores the diagonal elements (F _d ) of the flow matrix and memory 81b stores the diagonal elements (D _d ) of the distance matrix. With an asymmetric matrix, unlike with a symmetric matrix, these diagonal elements can contain non-zero values.

The register 82a holds one of f a, _a or f _b ,b read from the memory 81a, and at the timing when the other of f _a,a or f _b,b is read from the memory 81a, f _a,a or f One of _{b and b} is output and supplied to the difference calculation circuit 83a.

Register 82b holds either d φ _{(a), φ(a)} or d _{φ(b), φ(b)} read from memory 81b. Then, the register 82b stores d _{φ(a), φ(a)} or d _{φ (a), φ(a)} or d One _{of φ(b) and φ(b)} is output and supplied 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)} of the equation (44).

Multiplying circuit 84 calculates (f _b,b −f _a,a )(d _{φ(a), φ(a)} −d _{φ(b ), φ(b)} ).
The compensation term calculation circuit 85 calculates compensation terms (terms for compensation by eliminating the condition of skipping the calculation of i=a, b), which are the four items on the right side of equation (44). The compensation term calculation circuit 85 receives f _a,b , d _{φ(a), φ(b)} from the multiplexers 33a, 33b of the replica processing circuit 70a1 to calculate the compensation terms, and the multiplexer 33a of the replica processing circuit 70a2. , 33b receive f _b,a , d _{φ(b), φ(a)} .

Adder circuit 86 receives the value of one item on the right side of equation (44) from inner product calculation circuit 34 of replica processing circuit 70a1, and the value of two items on the right side of equation (44) from inner product calculation circuit 34 of replica processing circuit 70a2. receive. Further, adder circuit 86 receives the values of three items on the right side of equation (44) from multiplier circuit 84 and the values of four items on the right side of equation (44) from compensation term calculation circuit 85 . Then, the addition circuit 86 generates and outputs ΔC _asym by calculating the sum of these.

These circuits, memories, and the like can be included in the storage unit 11 or the processing unit 12 in FIG.
The configuration for storing the integer allocation vector φ representing the allocation state, the configuration for determining whether or not to accept the allocation change proposal that causes the calculated ΔC _ex , and the like are not shown in the drawing.

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

The above 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. As a recording medium, for example, a magnetic disk, an optical disk, a magneto-optical disk, a semiconductor memory, etc. can be used. Magnetic disks include flexible disks (FDs) and HDDs. Optical discs include CD (Compact Disc), CD-R (Recordable)/RW (Rewritable), DVD (Digital Versatile Disc) and DVD-R/RW. The program may be recorded on a portable recording medium and distributed. In that 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.
Computer 90 has CPU 91 , RAM 92 , HDD 93 , GPU 94 , input interface 95 , medium reader 96 and communication interface 97 . The 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 part of the programs and data stored in the HDD 93 into the RAM 92 and executes the programs. Note that the CPU 91 may include a plurality of processor cores in order to execute processing of a plurality of replicas in parallel, as shown in FIG. 4, for example. Computer 90 may also include multiple processors. A set of multiple processors (multiprocessor) may also be called a "processor".

The RAM 92 is a volatile semiconductor memory that temporarily stores programs executed by the CPU 91 and data used for calculation by the CPU 91 . The computer 90 may be provided with a type of memory other than the RAM 92, and may be provided with a plurality of memories.

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

The GPU 94 outputs an image (for example, an image representing a calculation result of an assignment problem, etc.) to a display 94a connected to the computer 90 according to instructions from the CPU 91 . As the display 94a, a CRT (Cathode Ray Tube) display, a liquid crystal display (LCD), a plasma display (PDP: Plasma Display Panel), an organic EL (OEL: Organic Electro-Luminescence) display, or the like can be used. .

The input interface 95 acquires an input signal from an input device 95 a connected to the computer 90 and outputs it to the CPU 91 . As the input device 95a, a mouse, a touch panel, a touch pad, a pointing device such as a trackball, a keyboard, a remote controller, a button switch, or the like can be used. Also, the computer 90 may be connected to a plurality of types of input devices.

The medium reader 96 is a reading device that reads programs and data recorded on the recording medium 96a. As the recording medium 96a, for example, a magnetic disk, an optical disk, a magneto-optical disk (MO), a semiconductor memory, or the like can be used. Magnetic disks include FDs and HDDs. Optical discs include CDs and DVDs.

The medium reader 96 copies the program and data read from the recording medium 96a to another recording medium such as the RAM 92 and the HDD 93, for example. The read program is executed by the CPU 91, for example. Note that the recording medium 96a may be a portable recording medium, and may be used for distribution of programs and data. Also, the recording medium 96a and the HDD 93 may be referred to as a computer-readable recording medium.

The communication interface 97 is an interface that is connected to a network 97a and communicates with other information processing apparatuses 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.

Although one aspect of the program, the data processing apparatus, and the data processing method of the present invention has been described above based on the embodiments, these are merely examples and are not limited to the above description.

For example, in the above description, the columns of the distance matrix are rearranged according to the allocation state. can get.

10 data processing device 11 storage unit 12 processing unit

Claims (9)

A program for causing a computer to execute a process of searching for a solution to an allocation problem by local search using an evaluation function representing a cost according to an allocation state,
out of the plurality of elements 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, stored in memory , calculating a first amount of change in the evaluation function when a first allocation change occurs in which the allocation destinations of the first element and the second element are exchanged, using vector arithmetic operations;
determining whether or not to allow the first allocation change based on the first amount of change;
When 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 exchanged. do,
A program that causes the computer to execute a process. If the assignment problem is QAP, allow the first assignment change based on a comparison result between the acceptance probability calculated based on the first change amount and the value of the temperature parameter and a random value. 2. The program according to claim 1, causing the computer to execute a process of determining whether or not. Before calculating the first amount of change,
when the assignment problem is QSAP, calculating a second amount of change in the evaluation function when a second assignment change occurs in which the first element is assigned to a first assignment destination;
determining whether or not to allow the second allocation change based on a comparison result between the acceptance probability calculated based on the second amount of change and the value of the temperature parameter and a random value;
updating the allocation state and the distance matrix if it is determined that the second allocation change is allowed;
After calculating the first amount of change,
Determining whether to allow the first allocation change based on a comparison result between the first amount of change and a predetermined value;
2. The program according to claim 1, which causes the computer to execute processing. reading the distance matrix row by row from the memory;
selecting two values of the two columns contained in the read row;
swapping storage locations of the two values in the memory;
By repeating the process, the two columns are exchanged,
2. The program according to claim 1, which causes the computer to execute processing. the memory further stores a transposed matrix of the distance matrix;
In the transposed matrix, among two rows corresponding to the two columns of the distance matrix, the first row is stored in a first shift register and the second row is stored in a second shift register;
The two columns are exchanged by repeating a process of swapping storage locations of two values output one by one from each of the first shift register and the second shift register and writing the values into the memory. I do,
2. The program according to claim 1, which causes the computer to execute processing. 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;
selecting two values of the two columns contained in the read row;
swapping storage locations of the two values and writing them to the second memory;
By repeating the process, the two columns are exchanged,
2. The program according to claim 1, which causes the computer to execute processing. performing the local search using parallel tempering with a plurality of replicas each having a different temperature parameter value;
the memory stores the distance matrix for each of a first replica and a second replica among the plurality of replicas;
calculating the first variation based on the distance matrix of the second replica while updating the distance matrix of the first replica;
2. The program according to claim 1, which causes the computer to execute processing. A data processing device that searches for a solution to an allocation problem by local search using an evaluation function representing a cost according to an allocation state,
a storage unit that stores 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;
A first evaluation function of the evaluation function when a first allocation change occurs in which the allocation destinations of the first element and the second element among the plurality of elements are exchanged based on the flow matrix and the distance matrix. using vector arithmetic operations, determining whether or not to permit the first allocation change based on the first amount of change, and determining to permit the first allocation change If so, a processing unit that 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 exchanged;
A data processing device having A data processing method for searching for a solution to an allocation problem by local search using an evaluation function representing a cost according to an allocation state,
the computer
out of the plurality of elements 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, stored in memory , calculating a first amount of change in the evaluation function when a first allocation change occurs in which the allocation destinations of the first element and the second element are exchanged, using vector arithmetic operations;
determining whether or not to allow the first allocation change based on the first amount of change;
When 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 exchanged. do,
Data processing method.

Images