Program, data processing device and data processing method

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

[Technical Field]

[0001] The present invention relates to a program, a data processing device, and a data processing method. [Background technology]

[0002] The assignment problem is a class of NP-hard combinatorial optimization problems with various 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).

[0003] 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 assigning n elements (facilities, etc.) to n destinations. In other words, the quadratic assignment problem is a problem of searching for an assignment that minimizes the value of the evaluation function expressed in the following equation (1). The evaluation function represents the cost according to the assignment state, and is also called the cost function.

[0004]

number

[0005] In equation (1), f i,j is the flow amount between elements with identification numbers i and j, d φ(i),φ(j) indicates the distance between the assignees to which the elements with identification numbers = i, j are assigned. By the way, there is a Boltzmann machine (also called an Ising machine) that uses an Ising evaluation function as a device for calculating large-scale discrete optimization problems, which von Neumann computers are not good at. The Boltzmann machine is a type of recurrent neural network.

[0006] A Boltzmann machine converts a combinatorial optimization problem into an Ising model that represents the spin behavior of a magnetic material. Then, the Boltzmann machine 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 the Boltzmann machine can also search for the state that maximizes the value of the evaluation function by changing the sign of the evaluation function. 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 0 or 1.

[0007] The Ising-type evaluation function is defined, for example, by the following equation (2).

[0008]

number

[0009] The first item on the right side is the sum of the values ​​(0 or 1) of two state variables selected from N state variables for all combinations of all state variables in the Ising model, without omissions or overlaps, and the weight value (which represents the strength of the interaction between the two state variables). i is the state variable with identification number i, s j is the state variable with identification number j, and w i,j is a weight value that indicates the magnitude of the interaction between the state variables with identification numbers i and j. The two terms on the right 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 the identification number = i.

[0010] Also, s i The change in energy due to the change in the value of i ) is expressed by the following equation (3).

[0011]

number

[0012] In equation (3), s i When changes from 1 to 0, Δs i becomes -1, and s i When changes from 0 to 1, Δs i is 1. Note that h i is called the local field, and Δs i Depending on h i Multiplying by the sign (+1 or -1) is ΔE i This becomes:

[0013] And, for example, ΔE i If is smaller than the noise value (also called thermal noise) obtained based on the values ​​of the random number and temperature parameters, then s i The process of updating the value of , which causes a state transition and updates the local field is repeated.

[0014] 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).

[0015]

number

[0016] In equation (4), x is a vector of state variables, and represents the allocation state of n elements to n allocation destinations. T is (x 1,1 ,…,x 1,n ,x 2,1 ,…,x 2,n ,……,x n,1 ,…,x n,n ) can be expressed as x i,j =1 indicates that the element with identification number = i is assigned to the element with identification number = j, and x i,j =0 indicates that the element with identification number=i is not assigned to the element with identification number=j.

[0017] W is a matrix of weight values, and is the flow quantity (f i,j ) and the matrix D of the distances between n assignees, it can be expressed by the following equation (5).

[0018]

number

[0019] [Patent Document 1] U.S. Patent Application Publication No. 2021 / 0326679 [Non-patent literature]

[0020] [Non-Patent Document 1] Eugene L. Lawler, “The quadratic assignment problem”, Management Science, Vol.9, No.4 pp.586-599, July 1963 [Non-patent document 2] 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 [Non-patent document 3] 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 [Non-patent document 4] 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 [Non-patent document 5] 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 [Non-patent document 6] 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 [Non-Patent Document 7] 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 [Non-patent document 8] 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 [Non-Patent Document 9] 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 [Non-Patent Document 10] 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 Summary of the Invention [Problem to be solved by the invention]

[0021] The method of calculating the quadratic assignment problem using the evaluation function described above requires a large number of repetitions of processes such as calculating the amount of change in energy and updating the local field, and also requires a large number of memory accesses, which makes the calculation time-consuming.

[0022] In one aspect, the present invention aims to provide a program, a data processing device, and a data processing method that can calculate an assignment problem at high speed. [Means for solving the problem]

[0023] In one embodiment, a program is provided that causes a computer to execute a process of searching for a solution to an allocation problem by local search using an evaluation function that represents a cost according to the allocation state, the program calculating a first change in the evaluation function using vector arithmetic operations when a first allocation change occurs in which the allocation 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 to be assigned to multiple destinations and a distance matrix that represents the distance between the multiple destinations, both stored in memory; determining whether or not to allow the first allocation change based on the first change; and if it is determined that the first allocation change is allowable, updating the allocation state and updating the distance matrix so that the two columns or two rows corresponding to the first element and the second element are swapped.

[0024] Also provided in one embodiment is a data processing device. Also, in one embodiment, a data processing method is provided. [Effects of the Invention]

[0025] In one aspect, the present invention allows for fast computation of the assignment problem. [Brief explanation of the drawings]

[0026] [Figure 1] FIG. 10 is a diagram illustrating an example of calculating a QAP. [Figure 2] FIG. 10 is a diagram illustrating an example of rearrangement of a distance matrix when calculating a QAP and an example of a data processing device. [Figure 3] FIG. 10 is a diagram illustrating an example of update calculation of a cache matrix. [Figure 4] FIG. 1 illustrates an example of a solver system that performs parallel tempering. [Figure 5] 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. [Figure 6] 10 is a flowchart showing the overall processing flow of local search by parallel tempering. [Figure 7] 10 is a flowchart illustrating an example of the flow of a replica initialization process in the case of a QAP. [Figure 8] 10 is a flowchart illustrating an example of the flow of a replica search process in the case of QAP. [Figure 9] FIG. 10 is a diagram illustrating an example of an algorithm for searching for a solution to a QSAP by local search using parallel tempering. [Figure 10] 10 is a flowchart showing the flow of an example of a replica initialization process in the case of QSAP. [Figure 11] 10 is a flowchart showing the flow of an example of replica search processing in the case of QSAP. [Figure 12] This figure 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. [Figure 13] 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. [Figure 14] FIG. 10 illustrates an example of load balancing. [Figure 15] FIG. 10 is a diagram showing the evaluation results of the degree of speedup of calculation processing by load distribution. [Figure 16] FIG. 10 is a diagram illustrating an example of a measurement algorithm. [Figure 17] This figure shows the relative speedup measured for the VΔC, SAM, and BM$ methods, and the memory hierarchy occupied depending on the problem size. [Figure 18] FIG. 10 is a diagram illustrating an example of a ΔC generation circuit. [Figure 19] FIG. 10 is a diagram illustrating a first example of a hardware configuration for performing column interchange. [Figure 20] FIG. 10 is a diagram illustrating an example of column interchange. [Figure 21] FIG. 10 is a diagram illustrating a second example of a hardware configuration for performing column interchange. [Figure 22]FIG. 10 is a diagram illustrating a first modified example of the second example of the hardware configuration for performing column interchange. [Figure 23] FIG. 10 is a diagram illustrating a second modification of the second example of the hardware configuration for performing column interchange. [Figure 24] FIG. 10 is a diagram illustrating a third example of a hardware configuration for performing column interchange. [Figure 25] FIG. 10 is a diagram illustrating a modification of the third example of the hardware configuration for performing column interchange. [Figure 26] FIG. 10 is a diagram illustrating another example of a ΔC generating circuit. [Figure 27] FIG. 10 is a diagram illustrating an example of a hardware configuration for performing processing on two replicas. [Figure 28] FIG. 10 is a diagram illustrating another example of a ΔC generating circuit. [Figure 29] FIG. 10 is a diagram illustrating an example of a replica processing circuit. [Figure 30] FIG. 10 is a diagram illustrating an example of a replica processing circuit used in calculating a QAP using an asymmetric matrix. [Figure 31] FIG. 1 is a diagram illustrating an example of hardware of a computer that is an example of a data processing device. DETAILED DESCRIPTION OF THE INVENTION

[0027] 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) as an example of an assignment problem by local search. QAP, QSAP, and local search will be described below.

[0028] (QAP) Figure 1 shows an example of how to calculate QAP. QAP is a problem that involves finding an allocation that minimizes the sum of the products of the flow rates between elements and the distances between the destinations to which each element is assigned, when allocating n elements (such as facilities) to n destinations.

[0029] FIG. 1 shows an example in which four facilities with identification numbers 1 to 4 are assigned to four allocation destinations (L1 to L4). The flow matrix representing the flow amount between n elements is expressed by the following equation (6).

[0030]

number

[0031] The flow matrix (F) is an n-row, n-column matrix. i,j is the flow amount in row i, column j, and represents the flow amount between elements with identification numbers = i, j. For example, the flow amount between facility with identification number = 1 and facility with identification number = 2 in Figure 1 is f 1,2 This can be expressed as:

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

[0033]

number

[0034] The distance matrix (D) is an n-by-n matrix. k,l is the distance between the k rows and l columns, and represents the distance between the assignees of identification numbers = k, l. For example, in Figure 1, the distance between the assignee (L1) of identification number = 1 and the assignee (L2) of identification number = 2 is d 1,2 This can be expressed as:

[0035] The QAP is calculated by searching for an assignment that minimizes the above 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 a set Φ n is an element of Φ n is the set of all permutations of the set N = {1, 2, 3, ..., n}. x is a binary variable contained in the binary state matrix X. i,j is expressed by the following equation (8).

[0036]

number

[0037] Figure 1 shows an example in which a facility with identification number 1 is assigned to a facility with identification number 2, a facility with identification number 2 is assigned to a facility with identification number 3, a facility with identification number 3 is assigned to a facility with identification number 4, and a facility with identification number 4 is assigned to a facility with identification number 1. The integer assignment vectors φ are φ(1)=2, φ(2)=3, φ(3)=4, and φ(4)=1. The binary state matrix X is expressed as x 1,2 , x 2,3 , x 3,4 , x 4,1 are each 1, and the other x i,j is set to 0.

[0038] In a QAP, either or both of the flow matrix and the distance matrix may be symmetric, or both of the flow matrix and the distance matrix may be asymmetric. In this embodiment, we mainly focus on QAPs using symmetric matrices (diagonal elements are 0 (biasless)). This is because such QAPs represent the majority of instances and simplify calculations. However, a QAP using a symmetric matrix can be directly converted into a QAP using an asymmetric matrix.

[0039] (QSAP) QSAP is a modified version of QAP. In QSAP, the number of elements is not equal to the number of destinations. For example, it is allowed to assign multiple elements to each destination. In QSAP, the distance matrix is ​​expressed by the following equation (9).

[0040]

number

[0041] 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.

[0042]

number

[0043] b i,k represents a fixed 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.

[0044]

number

[0045] The allocation state of n elements to m allocation destinations is expressed as an integer allocation vector ψ(ψ∈[1,m] n ), or the binary state matrix S. The binary variables contained in the binary state matrix S are s i,j is expressed by the following equation (12).

[0046]

number

[0047] (local search) In local search, a candidate solution is searched for within a nearby state that can be reached by changing the current state. One method of performing local search on 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), φ(b)) of two elements (represented by identification numbers = a, b) can be expressed by the following Equation (13).

[0048]

number

[0049] As shown in equation (13), the change (ΔC ex) is generated by the multiply-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).

[0050]

number

[0051] In equation (14), ΔB ex indicates the change in the constant cost for allocating elements to each destination, and is expressed by the following equation (15).

[0052]

number

[0053] 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 reassigning an element with identification number = a from its current allocation destination (identification number = ψ(a)) to an allocation destination with identification number = l is expressed by the following equations (16) and (17).

[0054]

number

[0055]

number

[0056] In local search, based on the change calculated as above, it is decided whether to accept a proposal for an allocation change that will cause a change in the value of the evaluation function by that amount. The decision to accept a proposal is based on predefined criteria, for example, a greedy method. When a greedy method is used, an allocation change that reduces the cost (value of the evaluation function (energy)) is accepted. When a greedy method is used, the proposal acceptance probability (PAR) is high at the beginning of the search, but tends to approach 0 later because the search gets stuck at a local minimum of the evaluation function and no further improvement is found.

[0057] 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, in order to add randomness to the change of allocation, a temperature parameter (T) is used to calculate the acceptance probability (P acc ) can be used.

[0058]

number

[0059] As T increases towards infinity, P acc and PAR increases, 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.

[0060] The calculation of ΔC for the QAP and 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 speeded up by using vector arithmetic operations as follows:

[0061] (Vectorization of ΔC calculation) Hereinafter, this method may be referred to as the VΔC method. exThe calculation of can be easily vectorized using inner product calculations unless the calculation of i=a,b is skipped.

[0062] The data processing device of this embodiment eliminates this condition and ex The calculation of is a product-sum operation loop for all n elements, and an additional product of the flow amount and the distance is added as a compensation term, as shown in the following equation (19), thereby performing compensation by eliminating the condition.

[0063]

number

[0064] 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 the elements with identification numbers = a and b are swapped.

[0065]

number

[0066]

number

[0067]

number

[0068] Δ in Equation (20) b a F represents the difference vector between row b and row a of the flow matrix, and F b,* denotes the bth row of the flow matrix, and F a,* denotes row a of the flow matrix. The elements of the vector ΔF are arranged in the order of the original flow matrix, while the elements of the vector ΔD are ordered to correspond to the current allocation state. This is done by transposing the binary state matrix X(XT (written as) and (D φ(a),* -D φ(b),* ) is shown mathematically by multiplying with D φ(a),* denotes the a row of the distance matrix, and D φ(b),* indicates the bth row of the distance matrix.

[0069] 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 multiple elements can be placed 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 the following equation (23).

[0070]

number

[0071] Using such a vector ΔD and the vector ΔF shown in equation (20), ΔC s ex can be calculated using the following equation (24).

[0072]

number

[0073] 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. In addition, ΔC due to rearrangement corresponding to equation (17) s rel The vectorized form of the calculation can be expressed as the following equation (25).

[0074]

number

[0075] 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ΔC method can be implemented on processors such as multiple CPUs (Central Processing Units) or multiple GPUs (Graphics Processing Units) by using SIMD (Single Instruction / Multiple Data).

[0076] However, although the rearrangement of the elements of ΔD performed by the VΔC method can be performed in a time proportional to the number of elements (n), it is not efficient because it is performed every time ΔC is calculated by swapping the allocation destinations of two elements. Also, it may require a large amount of calculation cost.

[0077] In order to minimize the calculation 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. In addition, hereinafter, the distance matrix rearranged in this way is referred to as the state-aligned D matrix, and when the allocation problem is QAP, it is referred to as the D X , if the assignment problem is QSAP, then D S It is written as follows.

[0078] 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 .

[0079] 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 a flash memory. The storage unit 11 may also include an electronic circuit such as a static random access memory (SRAM) register.

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

[0081] The processing unit 12 can be realized by a processor, which is hardware such as a CPU, a GPU, or a 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 an 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.

[0082] The processing unit 12 repeats the process of switching 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. The allocation state in which the smallest value among 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).

[0083] Figure 2 shows an example of rearrangement of the distance matrix when calculating the QAP. The allocation state at a certain time step = t and the state rearrangement D matrix corresponding to the allocation state are expressed as φ (t), D X(t) If a change of allocation to swap the allocation destinations of the elements with identification numbers 1 and 3 is proposed and accepted, then at time step t+1, the φ (t+1) In this case, the state alignment matrix D is set to D to correspond to the accepted allocation state. X(t) D with the first and third columns swapped X(t+1) is generated.

[0084] If the proposed reassignment is not accepted, φ (t) and D X(t) will not be updated. In the case of QAP, the state alignment D matrix (D X ) is the transpose matrix X of the binary state matrix X that represents the current allocation state, as shown in the following equation (26): T by the original distance matrix (D). An example of a hardware configuration for swapping the columns of the state alignment matrix D in hardware will be described later.

[0085]

number

[0086] Vector ΔD X can be calculated without rearranging the elements of the vector ΔD using the following equation (27), X By substituting into the above equation (22), ΔC ex can be calculated.

[0087]

number

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

[0089]

number

[0090] Vector ΔD due to allocation exchange and relocation S can be calculated without rearranging the elements of the above ΔD using the following equations (29) and (30), respectively, and by substituting them into the above equations (24) and (25), ΔC s ex and ΔC s rel can be calculated.

[0091]

number

[0092]

number

[0093] In the VΔC method described above, the elements of vector ΔD are rearranged by operations proportional to the number of elements = n each time ΔC is calculated. In contrast, in the above method (SAM method), D is rearranged to match the allocation state. X and D S The reordering of is only done if the proposed reassignment is accepted. This reduces the number of operations required per iteration on average, making the computation more efficient and therefore faster to solve the assignment problem.

[0094] (Boltzmann machine caching method (comparison example)) The following method can be considered for vectorizing the ΔC calculation. This method calculates the partial ΔC value (the local field (h i This method is sometimes referred to as the Boltzmann machine caching method (BM$ method).

[0095] For each bit in the binary state matrix (X or S), a cached local field is used. When a reassignment proposal is accepted, n 2 The matrix (hereinafter referred to as the cache matrix) of the n rows and n columns of the local field is updated in a time proportional to n, but the calculation of ΔC can be performed in a time that is independent of n.

[0096] In the case of QAP, the cache matrix (H) is calculated by Equation (31) at the start of the search process as follows: 3 is generated in a time proportional to

[0097]

number

[0098] Using the cached local field, we can generate the inner product ΔF ΔD using equation (32). Then, by substituting the generated ΔF ΔD into the above equation (22), we can obtain ΔC ex can be calculated.

[0099]

number

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

[0101]

number

[0102] 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 rows corresponding to zero elements of the vector ΔF.

[0103] In the case of QSAP, the cache matrix (H s ) is generated by equation (34) at the start of the search process.

[0104]

number

[0105] Using the cached local field, we can generate the inner product ΔF ΔD using equation (35). Then, by substituting the generated ΔF ΔD into the above equation (24), we can obtain ΔC s ex can be calculated.

[0106]

number

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

[0108]

number

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

[0110]

number

[0111] If the relocation reassignment proposal is accepted, the cache matrix is ​​updated according to equation (38).

[0112]

number

[0113] The BM$ method described above requires the storage of a cache matrix, but as the problem scale increases, the memory capacity required to store the cache matrix also increases. This makes it difficult to use a relatively small memory with high-speed readout capacity. The SAM method does not require the storage of such a cache matrix.

[0114] (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 explanation 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.

[0115] Each dedicated core stores the search instance data (D X , D S , H, H s ) has a dedicated cache to hold FIG. 4 is a diagram illustrating an example of a solver system that performs parallel tempering.

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

[0117] In the case of QAP, for each replica, the integer allocation vector φ, the cache matrix (H) (in the case of the BM$ method), and the state alignment matrix D (D X ) (in the case of the SAM method), the value of the evaluation function (C) is retained.

[0118] In the case of QSAP, for each replica, the integer allocation vector ψ and the cache matrix (H s ) (in the case of the BM$ method), the state alignment D matrix (D S ) (in the case of the SAM method), the value of the evaluation function (C) is retained.

[0119] 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 temperature parameters (T). min From T max Hereinafter, such a temperature parameter value that increases stepwise may be referred to as a temperature ladder.

[0120] The parallel tempering controller 21 controls the values ​​of adjacent temperature parameters (T k , T k+1 ) is set between the replicas, the value of the evaluation function (C k , C k+1 ) and T k , T k+1 Based on this, the exchange acceptance probability SAP is expressed by the following equation (39), T k , T k+1 Exchange the values ​​of

[0121]

number

[0122] 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, state-ordered D matrix, and cache matrix are stored in the memory unit 11 in the single-precision floating-point format, for example. This is to enable the use of fused multiply-and-accumulate instructions when performing inner product calculations or updating of cache matrices when executing the VΔC method or SAM method.

[0123] Fig. 5 is a diagram showing an example of an algorithm for searching for a solution to a QAP by local search using parallel tempering. In the example algorithm shown in Fig. 5, the number of replicas for which local search is performed in parallel is 32. In Fig. 5, "(1)" and "(31)" shown on lines 16 and 18 represent the above-mentioned formulas (1) and (31).

[0124] The algorithm begins by initializing the temperature parameter values ​​and replica states. Then the solver system enters an optimization loop, performing I local search iterations in parallel across all replicas.

[0125] After the first local search, a replica exchange process is performed, and the predefined BKS (Best-Known-Solution) cost (C BKS ) is found or the timeout limit is reached, the optimization loop resumes.

[0126] In the local search of QAP, a loop is executed for a preset number of iterations (I), two elements (facility in the example of Figure 5) are selected, and the ΔC ex The value of is calculated using the selected method (VΔC, BM$, or SAM).

[0127] Next, the calculated ΔC ex and the current temperature parameter value of each replica, P in Eq. (18) acc is calculated. Next, a random value between 0 and 1 is generated and P acc A Bernoulli trial is performed to determine whether the proposed swap of allocations between the two selected elements is acceptable. If the proposal is accepted, the state (allocation state) is updated.

[0128] 6 is a flowchart showing the overall processing flow of the local search by parallel tempering. In FIG. 6, the overall processing flow of the local search by parallel tempering of the algorithm shown in FIG.

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

[0130] Next, an optimization loop (steps S14 to S16) is performed while incrementing i one by one from i = 0 to i < M - 1. In the optimization loop, a replica search process is performed (step S15). In the replica search process, for each replica, I local searches are performed using the value of the set temperature parameter.

[0131] Thereafter, a replica exchange process (step S17) is performed, and it is determined whether to end the search (step S18). For example, as described above, when the pre-set BKS cost (C BKS ) is found or the time-out limit is reached, it is determined that the search has ended and the search ends. If it is determined not to end the search, the optimization loop from step S14 is restarted.

[0132] Note that when the search ends, the data processing device 10 may output the search result (for example, the minimum value, C min , φ min ). The search result may be displayed, for example, on a display device connected to the data processing device 10, or may be transmitted to an external device.

[0133] Note that FIG. 6 shows an example of performing the processing in order for each replica. For example, with a processor having 32 cores, the processing of steps S11, S12, and S15 for 32 replicas can be executed in parallel.

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

[0135] First, the integer assignment vector φ is initialized (step S20). In the process of step S20, for example, the assignment destination (placement destination of the facility) of the elements is determined randomly. Then, the initial cost (C) is calculated by the calculation of equation (1) (step S21).

[0136] After that, based on the initialized φ and C, the minimum value (C min and φ min ) is initialized (step S22). Also, corresponding to the initial value of φ, the state alignment D matrix (D X ) is initialized (step S23). Then, the process returns to the flowchart shown in FIG. 6.

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

[0138] First, two elements (identified by identification numbers = a, b) that are candidates for swapping the assignment destinations are selected (step S31). Then, the value of ΔC ex (a, b) when the assignment destinations of the two elements are swapped is calculated by equation (22) (step S32).

[0139] Next, using the calculated ΔC ex and the current temperature parameter values of each replica, P acc in equation (18) is calculated (step S33). Then, it is determined whether P acc is greater than a random value rand() in the range from 0 to 1 (step S34).

[0140] P acc If it is determined to be >rand(), D X Swap columns a and b (D X *,a ,D X *,b ←D X *,b ,D X *,a ) is performed (step S35). Furthermore, the allocation status is updated by switching 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)) is performed (step S37).

[0141] Then, C <C min It is determined whether or not C <C min If it is determined that min ←C, φ min ←φ) is performed (step S39).

[0142] After the process in step S39 or in step S34, acc If it is determined that the value is not rand(), or if the value is C <C min If it is determined that this is not the case, 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.

[0143] 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 a QSAP by local search using parallel tempering.

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

[0145] Next, taking the case of executing the SAM method as an example, 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 using a flowchart.

[0146] FIG. 10 is a flowchart showing an example of the flow of the replica initialization process in the case of QSAP. First, the integer assignment vector ψ is initialized (step S50). In the process of step S50, for example, the assignment destination (placement destination of the facility) of the elements is randomly determined. Then, the initial cost (C) is calculated by the calculation of formula (11) (step S51).

[0147] After that, based on the initialized ψ and C, the minimum value (C min and ψ min ) is initialized (step S52). Also, corresponding to the initial value of ψ, the state alignment D matrix (D S ) is initialized by formula (28) (step S53). Then, the process returns to the flowchart shown in FIG. 6.

[0148] FIG. 11 is a flowchart showing an example of the flow of the replica search process 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 incrementing i one by one from i = 1 to i < I.

[0149] First, an element with identification number = a and the assignment destination with 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 by formula (25) (step S62).

[0150] Next, using the calculated ΔC s rel and the current temperature parameter values of each replica, P acc of formula (18) is calculated (step S63). And Pacc It is determined whether or not is greater than a random number value rand() ranging from 0 to 1 (step S64).

[0151] P acc If it is determined to be >rand(), D s The values ​​of column a of the distance matrix D 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 ) is performed (step S67).

[0152] Then, C <C min It is determined whether or not C <C min If it is determined that min ←C, ψ min ←ψ) is performed (step S69).

[0153] After the process in step S69 or in step S64, acc If it is determined that the value is not rand(), or if the value is C in the process of step S68 <C min If it is determined that this is not the case, the process of step S70 is carried out.

[0154] In the process of step S70, two elements (represented by identification numbers a and b) are selected as candidates for replacement of allocation destinations. And, ΔC when the allocation destinations of the two elements are swapped s ex The value of is calculated by equation (24) (step S71).

[0155] Next, ΔC s ex It is then determined whether or not the value is <0 (step S72). Note that a predetermined value (fixed value) may be used instead of 0. ΔC ex s If it is determined that D is <0, S Swap columns a and b (D S *,a,D S *,b ←D S *,b ,D S *,a ) is performed (step S73). Furthermore, the allocation status is updated by switching the allocation destinations (ψ(a),ψ(b)←ψ(b),ψ(a)) (step S74), and the value of the evaluation function (C) is updated (C←C+ΔC s ex ) is performed (step S75).

[0156] Then, C <C min It is determined whether or not C <C min If it is determined that min ←C, ψ min ←ψ) is performed (step S77).

[0157] After the process of step S77 or in the process of step S72, ΔC s ex If it is determined that the value is not <0, or if the value is C <C min If it is determined that this is not the case, the process from step S61 is repeated until i = I. If i = I, the process returns to the process of the flowchart shown in FIG.

[0158] 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-type calculation 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 parallel tempering of 64 replicas.

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

[0160] When the SAM method is used, D X Despite the additional computational cost of updating, the calculation speed is improved compared to the VΔC method. Furthermore, when the BM$ method is used, the calculation speed is improved by 2.57 times in geometric mean compared to the VΔC method.

[0161] (Improved 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 vector-type arithmetic processing, SIMD built-in functions of AVX (Advanced Vector eXtensions)2 were used.

[0162] Using the same procedure as above, the same 10 QAP instances were run using vector-based computation. Figure 13 shows 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 Figure 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, SAM, and BM$ methods.

[0163] For scalar-type arithmetic processing, vector-type arithmetic processing is, on average, almost twice as fast as the VΔC method and more than three times as fast as the BM$ method and the SAM method. In the VΔC method, as described above, the rearrangement of the elements of ΔD performed in the VΔC method occupies a significant portion of the execution time, so the merit of the SIMD built-in function is the least.

[0164] (Dynamic load distribution) In parallel tempering, the PAR is different between the replica with the lowest temperature set and the replica with the highest temperature set. The larger the value of T (the higher the temperature), the higher the PAR. This may cause a large execution time gap between the threads performing replica processing in the update processing of the SAM method and the BM$ method.

[0165] To reduce this, the data processing device 10 tracks the time for each temperature iteration within the SLS function and uses that time to scale the number of iterations executed at each temperature.

[0166] FIG. 14 is a diagram showing an example of load distribution. In FIG. 14, the calculation time of ΔC before the exchange of temperature parameters and the update processing time for 32 replicas with a temperature ladder of T1 to T32 (T1 < T2 <... < T32) set are shown.

[0167] When load distribution is not performed, the update processing time becomes longer for replicas with larger temperature parameter values, and shorter for replicas with smaller temperature parameter values. Therefore, longer idle times occur for replicas with smaller temperature parameter values.

[0168] On the other hand, when load distribution is performed, for example, based on the execution time (thread execution time) of the replica with T1 set, the number of iterations of ΔC calculation in other replicas is scaled.

[0169] This allows thread execution times to be roughly equal among replicas, reducing idle time, as shown in FIG. To quantitatively demonstrate the effectiveness of such load balancing, we evaluated the degree to which load balancing accelerates calculation processing using the 10 instances mentioned above as the calculation targets. The results are shown below.

[0170] 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) is used whether load balancing is performed or not.

[0171] 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 the different PAR gaps between the extreme values ​​of the temperature parameter due to the characteristics of the instance families.

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

[0173] [Table 1]

[0174] Two types of speed-up are shown: incremental speed-up and cumulative speed-up. The 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).

[0175] (Benchmark results) For the VΔ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.

[0176] To measure the benchmark scores of the QAP solvers that implement the above three methods, instances in the QAP Library Reference (see Non-Patent Document 5) were used. In addition, 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 (without a known optimal solution) shown in Non-Patent Document 7. To measure the benchmark scores of the QAP solvers, we used a set of instances introduced in Non-Patent Document 8.

[0177] (QAP benchmark)

[0178] [Table 2]

[0179] Table 2 shows the TtS for each instance of the QAP using solvers that implement three vector-based computation methods (VΔC, SAM with load balancing, and BM$). 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 as the baseline (1.00). Table 2 also shows the TtS and 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 can solve some difficult instances with a 100% success rate within a 5-minute timeout window.

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

[0181] 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 difference in the results in Table 2 compared to Table 1 is due to the instances used.

[0182] 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. Also, these evaluation results are the results 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.

[0183] (QSAP benchmark) Compared to QAP, there has been little previous disclosure on QSAP solvers. The main reference is the PMITS (Parallel Memetic Iterative Tabu Search) algorithm implemented on a 20-core CPU (see Non-Patent Document 8). PMITS reaches the best solution with a timeout limit of 1 hour 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.

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

[0185] [Table 3]

[0186] As shown in Table 3, compared to QAP, the performance of SAM and BM$ methods compared to VΔC method across QSAP instances is almost 2 times worse, which can be attributed to their higher PAR.

[0187] (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.

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

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

[0190] [Table 4]

[0191] [Table 5]

[0192] Despite the short execution time, we found that we could improve BKS for all instances except four using SAM and BM$ methods, and the average cost was also lower than the previous BKS.

[0193] (Scaling analysis) Based on the benchmark results and a qualitative comparison between the VΔC and SAM methods, the SAM method appears to be more efficient. The VΔC method sequentially rearranges each element of ΔD at each iteration, while the SAM method rearranges D to match the allocation state. X and D S This is because the reordering of the rows will only occur if the proposed reassignment is accepted.

[0194] On the other hand, the relative performance of SAM and BM$ methods depends primarily on the 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.

[0195] Below are the results of a comparison of execution times depending on 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 ) is executed and the execution time of one loop is measured. At this time, a proposal with the desired PAR value is randomly accepted.

[0196] The problem sizes are divided into three groups based on the memory hierarchy used to store the replica data of the SAM method and the BM$ method. For problem sizes where n < 256 and 256 < n < 1024, I limit uses 100M and 10M respectively. For problem sizes where n > 1024, I limit uses 1M.

[0197] In 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 the performance while the CPU pipeline, cache, and memory are under full load, 64 loop instances were executed 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.

[0198] Since the density of the flow matrix affects the update function of the cache matrix (H) of the BM$ method, two separate simulation sets were executed. 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 by combining 10 problem sizes in the range of n = [100, 5000] and 29 PARs in the range of PAR = [0.0001, 0.1].

[0199] Figure 17 is a diagram showing the degree of relative speedup measured for the VΔC method, SAM method, and 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). Figure 17 shows the degree of speedup of the BM$ method with respect to the SAM method and the VΔC method for flow matrices of different densities. Furthermore, Figure 17 shows the degree of speedup of the SAM method with respect to the VΔC method. In the example of Figure 17, the memory hierarchy, in ascending order of storage capacity, is the L2 cache, L3 cache, and DRAM.

[0200] The maximum PAR values ​​for a given flow matrix density in a non-load-balanced parallel tempering simulation for 10 instances of QAP are shown in Table 6 below.

[0201] [Table 6]

[0202] The QSAP simulation was omitted because it is similar to the QAP simulation. There is no significant difference in the results of the two 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 (256kB in this example), and a group of four cores share an L3 cache (16MB in this example).

[0203] As shown in Figure 17, in the case of the SAM and BM$ methods, the D X The cache matrix (H) fits into the L2 cache of a core up to size n=256, and problems up to n=1024 fit into the L3 cache.

[0204] As shown in Figure 17, the relative performance (speedup) between the BM$ method and the other two methods decreases as the problem size increases and a memory hierarchy with a larger storage capacity is used. When the data used for search is moved 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.

[0205] For a solver that performs parallel tempering, as shown in Table 6, the maximum PAR value for the entire replica does not necessarily decrease with the problem size of the entire instance family. This indicates that for smaller QAP instances, the relationship between PAR and problem size to maintain the same relative speed as shown in Table 2 does not necessarily depend on the relationship as shown in Figure 17.

[0206] The degree of acceleration of the SAM method with respect to the VΔC method is divided into two groups based on the problem size. When the problem size is n ≤ 800 and most of the D for each loop instance fits in the cache, the SAM method is clearly more advantageous than VΔC in the range of PAR < 10%. Shifting the storage location of the data used for exploration from the L2 cache to the L3 cache has not significantly affected the relative performance between the two methods. X の大部分がキャッシュに収まる場合、SAM法では、PAR<10%の範囲ではVΔCよりも明らかに有利である。探索に用いるデータの保存先をL2キャッシュからL3キャッシュに移行することは、2つの方法間の相対的な性能に大きな影響を与えていない。

[0207] (Hardware Example) Hereinafter, a hardware example for implementing the SAM method will be described. For simplicity of explanation, hereinafter, it is assumed that both the flow matrix and the distance matrix are symmetric matrices (diagonal components are 0 (bias-free)). As described above, such QAPs are most of the instances, and this is for simplifying the calculation. A QAP using a symmetric matrix can be directly converted into a QAP using an asymmetric matrix.

[0208] The evaluation function of QAP using only symmetric matrices can be expressed as in the following equation (40).

[0209]

Equation

[0210] The difference between Equation (1) and Equation (40) is only that "j = 1" in Equation (1) becomes "j = i". When using the evaluation function represented by Equation (40), the ΔC used in the calculation of QAP excan be expressed by the following equations (41) and (42) instead of equations (19) and (22).

[0211]

number

[0212]

number

[0213] FIG. 18 is a diagram illustrating an example of a ΔC generating circuit. The ΔC generation circuit 30 includes a flow matrix memory 31a, a state-ordered D matrix memory 31b, difference calculation circuits 32a and 32b, multiplexers 33a and 33b, 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. 2.

[0214] 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 ) to remember.

[0215] 18, the flow matrix memory 31a and the state-ordered D matrix memory 31b are dual-port memories having two ports. Such dual-port memories can be used in some FPGAs.

[0216] The difference calculation circuit 32a calculates Δ b a F is F b,* -F a,* Calculate. The difference calculation circuit 32b calculates Δ φ(a) φ(b) D X D X φ(a),* -D X φ(b),* Calculate.

[0217] Multiplexer 33a is F a,* From f a,b Select and output. Multiplexer 33b is D X φ(a),* From d φ(a),b Select and output. The inner product calculation circuit 34 calculates Δ b a F and Δ φ(a) φ(b) D X The inner product calculation circuit 34 can be realized by, for example, a plurality of multipliers connected in parallel.

[0218] The register 35 holds the coefficient "2" included in the equation (42). Multiplication circuit 36 ​​is 2f a,b d φ(a),b Calculate. The adder circuit 37 adds 2f to the result of the inner product. a,b d φ(a),b By adding ΔC ex Calculate and output.

[0219] ΔC using such hardware ex The calculation is controlled by, for example, a control circuit (not shown) included in the processing unit 12 of FIG. 2 executing a program. As mentioned above, in the SAM method, ΔC ex When a swap of the allocation destinations of elements that generates (allocation change) is proposed and accepted, the columns of the state alignment D matrix are swapped to correspond to the accepted allocation states.

[0220] FIG. 19 is a diagram illustrating a first example of a hardware configuration for performing column interchange. 19, the columns of the state ordering D matrix stored in the state ordering 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.

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

[0222] The switch 41 writes the value output from the multiplexer 40a into the storage location in the state-ordered D-matrix memory 31b, where the value output from the multiplexer 40b was stored. The switch 41 also writes the value output from the multiplexer 40b into the storage location in the state-ordered D-matrix memory 31b, where the value output from the multiplexer 40a was stored.

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

[0224] First d 1,3 and d 1,1 are selected by multiplexers 40a and 40b, and the memory locations are swapped by switch 41. Next, d 2,3 and d 2,1 are selected by multiplexers 40a and 40b, and the memory locations are swapped by switch 41. The same process is repeated a total of n cycles to complete the column swapping.

[0225] Fig. 21 is a diagram showing a second example of a hardware configuration for performing column interchange. In Fig. 21, the same elements as those shown in Fig. 19 are denoted by the same reference numerals. As shown in FIG. 21, the state alignment D matrix memory 31b stores the state alignment D matrix as well as the transposed matrix ((D X ) T) are also stored. The columns of the state-ordered 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.

[0226] In the circuit configuration of FIG. 21, two rows of the transposed matrix corresponding to the two columns of the state sorting D matrix to be swapped are read out from the state sorting 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.

[0227] The shift register 45b holds n values ​​in the second 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.

[0228] The switch 41 switches the storage location in the state-ordered 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-ordered 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.

[0229] This hardware configuration also allows for column swapping in n cycles. Furthermore, the hardware configuration of Fig. 21 may have improved wiring efficiency compared to the hardware configuration of Fig. 19. Furthermore, the multiplexers 40a and 40b shown in Fig. 19, which may result in a relatively large circuit scale, are no longer necessary.

[0230] Fig. 22 is a diagram showing a first modification of the second example of the hardware configuration for performing column interchange. In Fig. 22, the same elements as those shown in Fig. 21 are denoted by the same reference numerals.

[0231] In the hardware configuration of FIG. 22, the transpose matrix ((D X ) T 21 in that the above two rows of the transposed matrix are stored in the transposed matrix memory 46. The other configurations are the same as those in FIG. 21.

[0232] Fig. 23 is a diagram showing a second modification of the second example of the hardware configuration for performing column interchange. In Fig. 23, the same elements as those shown in Fig. 22 are denoted by the same reference numerals.

[0233] In the hardware configuration of FIG. 23, in addition to the flow matrix, the flow matrix memory 31a shown in FIG. 18 stores the transposed matrix ((D X ) T ) are also stored. The columns of the state-ordered 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 out from the flow matrix memory 31a.

[0234] Fig. 24 is a diagram showing a third example of a hardware configuration for performing column interchange. In Fig. 24, 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.

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

[0236] The state alignment D matrix after the column interchange is copied to the state alignment D matrix memory 31b1. Even with this hardware configuration, it is possible to swap columns in n cycles. Also, the hardware configuration of Fig. 24 may have improved wiring efficiency compared to the hardware configuration of Fig. 19.

[0237] Fig. 25 is a diagram showing a modification of the third example of the hardware configuration for performing column interchange. In Fig. 25, the same elements as those shown in Fig. 19 are denoted by the same reference numerals. In the hardware configuration of FIG. 25, the flow matrix memory 31a shown in FIG. 18 stores the transposed matrix ((D X ) T ) is also stored. The columns of the state-ordered 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 out from the flow matrix memory 31a.

[0238] 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. 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. a,* and stores the flow matrix memory 51a to F b,* When F is read a,* A register 52a is used to output the difference signal D and supply it to the difference calculation circuit 32a and the multiplexer 33a. X φ(a),* and stores the state-aligned D matrix memory 51b to D X φ(b),* When D is read X φ(a),* and supplies it to the difference calculation circuit 32b and multiplexer 33b.

[0239] 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.

[0240] Fig. 27 is a diagram illustrating an example of a hardware configuration for processing two replicas, which combines the ΔC generation circuit 50 shown in Fig. 26 with the hardware configuration for column exchange shown in Fig. 22.

[0241] In the example of FIG. 27, the state alignment D matrix memory 51b stores the state alignment D matrix (D X R1 , D X R2 ) is stored. The assignment change is accepted by one replica, and a column exchange is performed for the one replica using the write port of the state-ordering D-matrix memory 51b. Meanwhile, for the other replica, ΔC ex The process of calculating is performed.

[0242] If neither replica is performing an update, the process for one replica is performed by subtracting ΔC ex is stalled while the calculation of

[0243] 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. Only symmetric permutation problems are calculated, and ΔC ex is expressed by the following equation (43), a ΔC generation circuit 60 as shown in FIG. 28 can also be used.

[0244]

number

[0245] 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 calculates d φ(a),i -d φ(b),i The selector 61 has n two-input, one-output multiplexers 61ai that select and output either the output of the subtractor 32bi or 0.

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

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

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

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

[0250] In addition, the configuration for storing the integer allocation vector φ representing the allocation state and the calculated ΔC ex The configuration for determining whether or not to accept a proposal for an allocation change that would cause this is not shown in the figure.

[0251] (Extension to QAP when using asymmetric matrices) In the calculation of QAP using an asymmetric matrix (diagonal elements are non-zero), ΔC ex ΔC equivalent to asym is expressed by the following equation (44).

[0252]

number

[0253] Such a ΔC asym The replica processing circuit that calculates ΔC ex This can be realized using a replica processing circuit that calculates Fig. 30 is a diagram showing an example of a replica processing circuit used in calculating a QAP using an asymmetric matrix. In Fig. 30, the same elements as those shown in Fig. 29 are denoted by the same reference numerals.

[0254] ΔC asym The replica processing circuit 80 that calculates ΔC ex However, the flow matrix memory 51a of one of the replica processing circuits 70a1 and 70a2 stores the transposed matrix of the flow matrix (F T ) is stored in the state ordering D matrix memory 51b of the replica processing circuit 70a2, and the transposed matrix ((D X ) T ) is stored.

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

[0256] The memory 81a stores the diagonal elements of the flow matrix (F d ) and the memory 81b stores the diagonal elements of the distance matrix (D d ) When using a non-symmetric matrix, unlike when using a symmetric matrix, these diagonal elements can contain non-zero values.

[0257] The register 82a stores the f a,a or f b,b and stores one of the following in memory 81a to f. a,a or f b,b When the other data is read, a,a or f b,b One of the two is output and supplied to the difference calculation circuit 83a.

[0258] Register 82b is read from memory 81b. φ(a),φ(a) or d φ(b),φ(b) The register 82b holds one of the two values ​​from the memory 81b. φ(a),φ(a) or d φ(b),φ(b) When the other data is read, φ(a),φ(a) or d φ(b),φ(b) One of the two is output and supplied to the difference calculation circuit 83b.

[0259] The difference calculation circuit 83a calculates f in equation (44). b,b -f a,a Calculate. The difference calculation circuit 83b calculates d in equation (44). φ(a),φ(a) -d φ(b),φ(b) Calculate.

[0260] The multiplication circuit 84 calculates (f b,b -f a,a )(d φ(a),φ(a) -d φ(b),φ(b) ) is calculated. The compensation term calculation circuit 85 calculates the compensation terms (the four items on the right side of equation (44) that compensate for the elimination of the condition that the calculation of i=a, b is skipped). In order to calculate the compensation terms, the compensation term calculation circuit 85 uses multiplexers 33a and 33b of replica processing circuit 70a1 to calculate the compensation terms. a,b , d φ(a),φ(b) and outputs f from the multiplexers 33a and 33b of the replica processing circuit 70a2. b,a , d φ(b),φ(a) Receive.

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

[0262] These circuits, memories, etc. may be included in the storage unit 11 or the processing unit 12 in FIG. In addition, the configuration for storing the integer allocation vector φ representing the allocation state, the calculated ΔC ex The configuration for determining whether or not to accept a proposal for an allocation change that would cause this is not shown in the figure.

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

[0264] 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.

[0265] 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.

[0266] 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, as shown in FIG. 4, for example. The computer 90 may also include multiple processors. A collection of multiple processors (multiprocessor) may also be called a "processor."

[0267] The RAM 92 is a volatile semiconductor memory that temporarily stores programs executed by the CPU 91 and data used in calculations by the CPU 91. Note that the computer 90 may be provided with a type of memory other than the RAM 92, or may be provided with multiple memories.

[0268] The HDD 93 is a non-volatile storage device that stores software programs such as an 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.

[0269] The GPU 94 outputs an image (for example, an image showing the calculation result of an assignment problem) to a display 94a connected to the computer 90 in accordance with an instruction from the CPU 91. The display 94a may be a CRT (Cathode Ray Tube) display, a liquid crystal display (LCD: Liquid Crystal Display), a plasma display (PDP: Plasma Display Panel), an organic EL (OEL: Organic Electro-Luminescence) display, or the like.

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

[0271] The medium reader 96 is a reading device that reads programs and data recorded on a recording medium 96a. The recording medium 96a may be, for example, a magnetic disk, an optical disk, a magneto-optical disk (MO: Magneto-Optical disk), or a semiconductor memory. Magnetic disks include FDs and HDDs. Optical disks include CDs and DVDs.

[0272] The medium reader 96 copies programs and data read from the recording medium 96a to another recording medium such as the RAM 92 or the 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 the HDD 93 may also be referred to as computer-readable recording media.

[0273] The communication interface 97 is an interface connected to a 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 may be a wireless communication interface connected to a base station via a wireless link.

[0274] While one aspect of the program, data processing device, and data processing method of the present invention has been described above based on the embodiment, these are merely examples and the present invention is not limited to the above description.

[0275] For example, in the above explanation, the columns of the distance matrix are rearranged according to the allocation state, but the same effect can be obtained by rearranging the rows of the distance matrix according to the allocation state by modifying the formula appropriately. [Explanation of symbols]

[0276] 10 Data processing device 11 Storage section 12 Processing section

Claims

1. 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.

2. 2. The program according to claim 1, which causes the computer to execute a process of determining whether to accept the first allocation 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 when the allocation problem is QAP.

3. 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.

4. 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.

5. 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.

6. 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.

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

8. 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:

9. 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.

Images