Koopman analysis of combinatorial optimization problems with replica exchange Monte Carlo method
Tatsuya Naoi Tatsuya Kishimoto Jun Ohkubo






Graduate School of Science and Engineering, Saitama University, Sakura, Saitama 338–8570 Japan

§ ABSTRACT 
Combinatorial optimization problems play crucial roles in real-world applications, and many studies from a physics perspective have contributed to specialized hardware for high-speed computation. However, some combinatorial optimization problems are easy to solve, and others are not. Hence, the qualification of the difficulty in problem-solving will be beneficial. In this paper, we employ the Koopman analysis for multiple time-series data from the replica exchange Monte Carlo method. After proposing a quantity that aggregates the information of the multiple time-series data, we performed numerical experiments. The results indicate a negative correlation between the proposed quantity and the ability of the solution search.


Koopman analysis of combinatorial optimization problems with replica exchange Monte Carlo method
    Tatsuya Naoi, Tatsuya Kishimoto, and Jun Ohkubo
    September 9, 2024
================================================================================================





§ INTRODUCTION



Combinatorial optimization problems play a crucial role in various fields in real-world applications, and their solutions contribute to efficient resource utilization, cost reduction, and service quality improvement. Examples include supply chain management, finance, manufacturing, transportation, energy, healthcare, and telecommunications. In combinatorial optimization problems, we find combinations of variables that minimize an objective function while satisfying various constraints. Since the minimization procedure is the same as seeking ground energy states in statistical physics <cit.>, there are many works for combinatorial optimization problems from a physics viewpoint. For example, D-wave Advantege by D-Wave systems <cit.> employs the principle of quantum annealing <cit.>. Another type of hardware is based on complementary metal-oxide semiconductors (CMOS), including devices developed by Hitachi <cit.>, Toshiba <cit.>, and Fujitsu <cit.>. In the CMOS-type hardware, an architecture based on the replica exchange Monte Carlo method <cit.> is sometimes employed, which exploits the merit of parallel computation.


While recent annealing hardware enables us to solve combinatorial optimization problems at high speed, the time required to find the optimal solution varies from problem to problem. For example, a class of 0-1 quadratic knapsack problem <cit.> is NP-hard, but some are easy to solve, and others are not. Conventional simulated or quantum annealing requires a slow reduction in temperature or quantum effects when the combinatorial optimization problem is hard to solve. Although there is a well-known schedule that yields an optimal solution <cit.>, the schedule is too late in practice. In addition, it is worth knowing in advance whether the problem is easy to solve or not.

Here, we expect that the time-series data used to find the solution will reflect the degree of difficulty in problem-solving; the dynamics for easy problems will have simple structures, while difficult problems yield complicated dynamics. To investigate the dynamics, one can employ the Koopman theory <cit.>, which has attracted attention in various fields. The Koopman theory enables us to use linear algebra for nonlinear systems; it is also possible to analyze stochastic nonlinear systems within the framework of linear algebra. In the Koopman theory, we deal with a Koopman operator for the time evolution in the observable function space, not the state space. While the observable function space is infinite-dimensional, the Koopman operator is linear even if the underlying dynamical system is nonlinear. Several methods have been developed in recent years to obtain Koopman operators from an observed data set. One of the methods is dynamic mode decomposition (DMD), which analyzes dynamical systems from observed data <cit.>. There are extended methods, such as extended dynamic mode decomposition (EDMD) <cit.>. The EDMD can efficiently approximate the Koopman operator as a finite-dimensional matrix. Applications of EDMD are currently being studied, including time series prediction, system identification <cit.>, and control <cit.>. 


In this paper, we aim to clarify the relationship between the Koopman analysis and the degree of difficulty in combinatorial optimization problems. The EDMD yields an approximated matrix, the so-called Koopman matrix, for the Koopman operator. Since the replica exchange Monte Carlo method generates time-series data for each temperature, we get several datasets for the Koopman analysis. After applying the EDMD to time-series data generated from the replica exchange Monte Carlo method, we analyze the eigenvalues of the derived Koopman matrix for each temperature; the information of the eigenvalues will contain that of the degree of difficulty for each problem. We will also discuss a conventional data analysis based on singular value decomposition (SVD).


This paper is organized as follows. Section <ref> reviews the background of the annealing method and the Koopman theory. Section <ref> yields the main proposal, including data acquisition and analysis methods. In Sect. <ref>, we demonstrate the proposed method with numerical experiments on randomly generated problems and the 0-1 quadratic knapsack problems; the comparison with the proposed method and the conventional SVD are also discussed. We give a summary and mention of future work in Sect. <ref>.





§ BACKGROUND





 §.§ QUBO formulation and Ising model


A quadratic form of binary variables called quadratic unconstrained binary optimization (QUBO) formulation is widely used as an input for specialized hardware. The QUBO formulation is deeply related to the Ising model, and both have binary variables. The QUBO formulation has the following cost function for the state vector a:

    E(a) = ∑_i∈𝒟∑_j∈𝒟1/2 Q_ija_ia_j,

where 𝒟 is the set of indices of the variables, a_i∈{0,1} is the i-th binary variable in a, and Q_ij∈ℝ is the strength of the interaction between the binary variables a_i and a_j. The number of spins is N, i.e.,  |𝒟| = N. By contrast, the Ising model has binary spins with σ_i∈{-1,1}, and the energy (cost) function of the Ising model is denoted as follows:

    E(σ) = -∑_i∈𝒟∑_j∈𝒟1/2 J_ijσ_iσ_j-∑_i∈𝒟h_iσ_i,

where σ is the spin vector, J_ij∈ℝ is the two-body interaction between the spins σ_i and σ_j, and h_i is the external magnetic field on σ_i.

The energy function of the Ising model is equivalent to the cost function of QUBO formulation via the following transformation:

    a_i = 1+σ_i/2.

Conversion between {J_ij},{h_i} and {Q_ij} is also possible.

It is possible to convert various combinatorial optimization problems into the QUBO formulation <cit.>. Hence, minimizing the functions in Eqs. (<ref>) or (<ref>) corresponds to seeking solutions to combinatorial optimization problems. 







 §.§ Replica exchange Monte Carlo method



The replica exchange Monte Carlo method, also known as the parallel tempering method, improves sampling efficiency in Monte Carlo simulations and Markov chain Monte Calro methods <cit.>. In the replica exchange method, the temperature of each replica is determined and fixed in advance, and each replica develops with the conventional Metropolis rule. There is an exchange procedure for replicas. Due to the replica exchanges, when a low-temperature replica gets stuck in a local minimum, the replica can escape from a local minimum by exchanging a higher-temperature replica. 

Let R be the number of replicas and a^(r) be the state vector for the r-th replica. Then, the Gibbs distribution is defined as follows:

    P(a^(r))=1/Z^(r)exp(-E(a^(r))/T^(r)),

where T^(r), E(a^(r)), and Z^(r) correspond to the temperature parameter, the cost function, and the normalization constant for the r-th replica, respectively. While it is difficult to evaluate the normalization constant Z^(r) in general, there is no need to evaluate it; the state transition obeys the following conventional Metropolis rule:

    P^(r)_change=min[1,exp(-Δ E^(r)/T_r)],

where Δ E^(r) is the energy difference from the previous state to the next one in the r-th replica. The Metropolis rule yields a sampling from the Gibbs distribution without evaluating the normalization constant Z^(r).

As denoted above, there is the exchange procedure of the state vectors between different replicas. Although low-temperature settings are necessary when seeking stable states, there are many local minima that have difficulty escaping at low temperatures. Hence, we seek various configurations using the replicas with high temperatures. Here, we define the temperature parameters in ascending order: T_1 < T_2 < ⋯ < T_R. Then, the probability P^(m,l)_exchange for the exchange between the m-th and l-th replicas is defined as

    P^(m,l)_exchange=min[1,exp{(E^(m)-E^(l))(1/T_m-1/T_l)}],

where E^(m) and E^(l) are the energies of the m-th and l-th replicas, respectively. In this paper, we consider only exchanges at two adjacent replicas, and the following equation holds in Eq. (<ref>):

    l = m + 1

for m=1, ⋯, R-1.

Note that the data analysis method proposed in Sect. <ref> assumes the usage of the replica-exchange Monte Carlo method. As denoted in Sect. <ref>, the replica-exchange Monte Carlo method is employed on a certain type of annealing hardware and is also suitable for parallel computing. Hence, the assumption of using the replica-exchange Monte Carlo method is appropriate.





 §.§ Koopman theory


We here briefly review the Koopman theory; for more details, see, for example, the review paper in Ref. <cit.>.

The Koopman operator is a linear operator defined mainly in nonlinear dynamical systems. In the Koopman theory, we consider the time evolution of an observable function instead of the state variables. Consider the following time-evolution:

    a_t+1 = F(a_t),

where a_t is a state vector at time t, and F is a nonlinear time evolution operator. Here, we consider a deterministic dynamics for simplicity; a stochastic case will be commented later. We also introduce an observable function ϕ(a). For example, ϕ(a) = a_i corresponds to the observation of the i-th spin variable state. Then, the Koopman operator 𝒦 acts on the function ϕ as follows:

    𝒦ϕ = ϕ∘ F.

Equations (<ref>) and (<ref>) lead to

    ϕ(a_t+1)=(𝒦ϕ)(a_t).

It is easy to confirm the linearity of the Koopman operator; for any constants c_1, c_2 ∈ℝ and observable functions ϕ_1, ϕ_2, we have

    {𝒦(c_1 ϕ_1 + c_2 ϕ_2) } (a_t)    = (c_1 ϕ_1+c_2 ϕ_2)(a_t+1) 
       = c_1 ϕ_1(a_t+1) + c_2 ϕ_2(a_t+1) 
       = c_1(𝒦ϕ_1)(a_t) + c_2(𝒦ϕ_2)(a_t) 
       = (c_1𝒦ϕ_1 + c_2𝒦ϕ_2)(a_t).


Note that the Koopman operator 𝒦 is infinite-dimensional because 𝒦 acts on elements in a function space.  Therefore, it is necessary to approximate the Koopman operator 𝒦 as a finite-dimensional Koopman matrix K in practice. To describe the Koopman matrix K, we introduce a so-called dictionary. The dictionary is defined as

    ψ(a) = [ψ_1(a),ψ_2(a),⋯,ψ_N_dic(a)]^⊤,

where ψ_i(a) is the i-th function and N_dic is the size of the dictionary. Note that in Ref. <cit.>, the dictionary is defined as a row vector, but in this paper we consider it as a column vector. Then, we introduce the following linear combination to express an observable function:

    ϕ(a_t) = ∑_k=1^N_dicc_kψ_k(a_t) = c^⊤ψ(a_t),

where c is a coefficient vector. Combining Eqs. (<ref>) and (<ref>), we have

    (𝒦ϕ)(a_t)  = c^⊤(𝒦ψ)(a_t).

Note that c is time-independent. Hence, instead of the the action of the Koopman operator on the observable function ϕ, it is enough to consider the action on the dictionary as follows:

    ψ(a_t+1) = 𝒦ψ(a_t) ≃ K ψ(a_t),

which leads to the Koopman matrix K.

The EDMD is a method to derive the Koopman matrix K from a dataset <cit.>. Here, we consider a single time-series dataset {a_1,a_2,…,a_M}, although it is sufficient to have snapshot pairs rather than the single time-series data. Then, the least-squares problem with the cost function,

    J = ∑_t=1^M ψ(a_t+1) - Kψ(a_t) ^2,

leads to the Koopman matrix K immediately. The solution is denoted as

    K = QG^+,

where 

    G = 1/m∑_t=1^Mψ(a_t)^⊤ψ(a_t),   
    Q = 1/m∑_t=1^Mψ(a_t+1)^⊤ψ(a_t),

and G^+ is the pseudo-inverse of the matrix G.

In this paper, the Metropolis rule in Eq. (<ref>) yields the dynamics, and hence we should consider a stochastic dynamics. As discussed in Ref. <cit.>, the Koopman matrix K gives expectations in the stochastic case; Eq. (<ref>) is replaced simply with

    𝔼[ψ(a_t+1)] ≃ K ψ(a_t).

Even in the stochastic case, the Koopman matrix K contains the information of the system dynamics.




§ DIFFICULTY IN PROBLEM-SOLVING AND KOOPMAN ANALYSIS



Here, we will discuss a relationship between the Koopman analysis and the degree of difficulty in problem-solving. 



 §.§ Flow for analyzing data


One would expect that the dynamics for easy problems will have simple structures, while difficult problems yield complicated dynamics. The Koopman matrix contains the degree of complexity of the dynamics; it is possible to describe the dynamics of simple structures with a small number of modes, and the dynamics of complex structures will require many modes. Therefore, the distribution of eigenvalues should reflect the information of the dynamics.

However, a single time series of data alone cannot determine whether the eigenvalue distribution is simple or not. For example, the state of any problem will hardly change at considerably low temperatures, and only a few modes will remain. Therefore, we focus on the variation of the eigenvalue distributions at different temperatures. If the difficulty in problem-solving is high, the number of modes needed to describe the dynamics will increase as soon as the temperature becomes high. 





From the above considerations, we employ several time-series datasets generated from the replica exchange Monte Carlo method. Figure <ref> shows the flow for analyzing data. The replica exchange Monte Carlo method generates several time-series data with different temperature parameters, and we simply analyze each of those datasets. Although the flow is simple, there are some notices for the data analysis. We explain each of those notices below.




 §.§ Reduction of dictionary size


The dictionary size in the EDMD becomes enormously large when the number of variables increases. For simplicity, consider monomial dictionary functions with the form a_1^p_1 a_2^p_2⋯ a_N^p_N where p_i ∈ℕ_0. When L-1 is the maximum degree of each variable, i.e., p_i < L for all i, the dictionary size is L^N. Of course, we can restrict L to be less than two because the spin only takes 0 or 1; for example, a_1^2 = a_1. However, we still have the dictionary size with 𝒪(2^N). Due to the curse of dimensionality, we cannot evaluate the Koopman matrix.

In the data analysis, we restrict the total degree of monomials to be less than two and remove a_i^2 for i = 1, …,N from the dictionary. Hence, we employ the following dictionary:

    ψ(a) =
    [
    1, a_1, a_2, …, a_N, a_1 a_2, a_1 a_3, …, a_N-1 a_N]^⊤.

Preliminary numerical verification confirms that a higher-order dictionary, which also includes monomials of degree 3 or higher, does not significantly affect the results.  We consider the reason for this is that the QUBO formulation has only up to two-body interactions. Hence, we employ the dictionary in Eq. (<ref>) in the data analysis.




 §.§ Preparation of snapshot pairs


As denoted in Fig. <ref>, there are replica exchanges between different temperature settings. In the EDMD, the analysis is based on snapshot pairs. Hence, one might suspect that the snapshot pairs at replica exchange timings could cause some problems in the data analysis; the snapshot pairs at the timings do not reflect the true dynamics of the system. 

For the preparation of the snapshot pairs, we performed preliminary numerical experiments in which we removed the snapshot pairs at the replica exchange timings. The numerical results showed no significant differences due to data preprocessing. The reason is the small number of replica exchanges; these snapshot pairs do not affect the data analysis. Hence, in the flow in Fig. <ref>, we do not remove the snapshot pairs at the replica exchange timings for simplicity.





 §.§ Summation of eigenvalues


After the flow in Fig. <ref>, we have several eigenvalue distributions. As denoted above, we investigate the variation of the eigenvalue distribution at different temperatures. 





Figure <ref>(a) shows examples of the eigenvalue distributions; the experimental setting is the same with Sect. <ref>; we will explain it later. In Fig. <ref>(a), the distribution varies depending on the problem. There are several ways to capture the characteristics of the eigenvalue distributions. For example, the number of eigenvalues above a certain threshold would represent the number of dominant modes. After trying some ways, we finally decided to employ the summation of the eigenvalues. Figure <ref>(b) plots the sum of eigenvalues versus each temperature; we see clear differences depending on the problem.




 §.§ Characterization of temperature dependency



Finally, we characterize the shape of the temperature dependency in Fig. <ref>(b). After some trials, we found that the area of the rapidly changing region of the temperature dependency in Fig. <ref>(b) captures the characteristics of the difficulty of problem-solving. For example, the left side of the dotted line in Fig. <ref>(b) corresponds to the region showing sharp temperature dependency for each case.

Let w_r be the sum of Koopman eigenvalues at the r-th temperature. We introduce the following quantity S to characterize the difficulty in problem-solving:

    S = ∑^R'_r=1(1/2(w_r+1+w_r)(T_r+1-T_r)),

where T_r is the temperature of the r-th replica, and R' is the index of the end temperature of the area. Although one can determine the index R' by the figure appearance, we here employ the following determination method:

  
  * Set r = 3.
  
  * Calculate Δ_r = (w_r+1-w_r)/(w_r-w_r-1).
  
  * If Δ_r < 0.5, then R'=r. If not, r → r+1 and go back to Step 2.

This procedure does not use the first few indices because the data does not show a monotonic increase in the first three or so temperature ranges. Additionally, this rule is not universally applicable; increasing the number of replicas and sampling more finely across temperature ranges could not necessarily exhibit a monotonic increase. However, in the numerical experiments below, we determine R' in the manner described above to automate the procedure.






§ NUMERICAL EXPERIMENT



We perform several numerical experiments to check whether the proposed method successfully captures the characteristics of the difficulty of problem-solving. In the following numerical experiments, we deal with two examples; one is randomly generated, and the other is the 0-1 quadratic knapsack problem <cit.>, which was also used in the previous work for the specialized hardware <cit.>. As for the 0-1 knapsack cases, we generated some problems randomly to vary the difficulty level of problem-solving.

The purpose of the following numerical experiments is to examine the relationship between the Koopman analysis and difficulty in problem-solving. Hence, we compare the number of finding times of the optimal solution with the quantity proposed in Sect. <ref>. Note that we consider only problems with small spin numbers because we need optimal solutions for the comparison.




 §.§ Randomly generated problems



Although we tried several problem settings, the following method for generating problems of varying difficulty is employed here. Here, the number of spins is 15, and we randomly generate J_ij in Eq. (<ref>) from 𝒩(3,0.5) or 𝒩(-3,0.5) with equal probability. The coefficients h_i in Eq. (<ref>) is also generated from 𝒩(2,0.5) or 𝒩(-2,0.5) with equal probability. The number of replica is R=15, and the temperatures {T_r} for r=1, 2, ⋯, R are defined as follows:

    T_r = r + 1.0.

We obtain 10,000 time-series data for each temperature by using the replica exchange Monte Carlo method. Before the analysis by the EDMD, the binary variables {-1,1} are converted to {0,1} by Eq. (<ref>). Here, we generate nine different problems, and for each problem, the quantity S by the data analysis flow described in Sect. <ref> is evaluated.






Figure <ref>(a) shows the temperature-dependency of the the sum of Koopman eigenvalues. For the evaluation of the quantity S in Eq. (<ref>), R' is determined by the procedure described in Sect. <ref> as follows: 

  
  * R'=5 for Cases 1, 2, 3, 4 and 7,
  
  * R'=6 for Cases 5 and 6,
  
  * R'=7 for Case 9,
  
  * R'=8 for Case 8.


Next, we measure the number of finding times of the optimal solution for each case. For the 15 spins problem, the number of states is 2^15 = 32,768, and hence, it is possible to search exhaustively to find the state that minimizes the cost function. Note that the temperature parameters used to solve the problem in the replica exchange Monte Carlo method differ from those used to compute the features S; we employ the following conventional settings:

    T_r = 0.001 + (r/R)^2,    (R=5  and  r=1, 2, ⋯, R).

We apply the updates for a randomly chosen spin N times and the replica exchange procedure at the timing; the replica exchange occurs between a randomly selected replica number and the number one above it. We call this procedure an iteration step. After repeating the iteration steps 100 times, we judge whether we find the optimal solution. We repeat the search 100 times and count the number of finding times of the optimal solution. Figure <ref>(b) shows the relationship between the number of finding times of the optimal solution in 100 trials and the feature S; it shows the results of ten runs of the above experiment for each problem case. 


In Fig. <ref>(b), we see a clear negative correlation between the number of finding times the optimal solution and the feature S. This numerical results indicate that our conjecture is valid: the temperature-dependency in the sum of Koopman eigenvalues captures the difficulty in problem-solving.





 §.§ Quadratic knapsack Problem



The 0-1 quadratic knapsack problem extends the standard knapsack problem. The standard knapsack problem aims to select items within a given capacity constraint to maximize the total value. However, in the 0-1 quadratic knapsack problem, the value of the items is represented in a quadratic form. As denoted above, the 0-1 quadratic knapsack problem was discussed in the specialized hardware previously <cit.>. The cost function is denoted as follows:

    E(x) = -1/2∑_i∈ D∑_j∈ DVal_ijx_ix_j - ∑_i∈ DVal_ix_i
    
    +λmax(0,∑_i∈ DWT_ix_i-WT_max),

where Val_i and WT_i mean the value and the weight of the i-th item, respectively; Val_ij is the interaction value between the i-th item and j-th item, and WT_max is the capacity of the knapsack. Note that x_i∈{0,1}. 




Here, we consider problems with 15 items. The problem parameters,  Val_ij, Val_i, WT_i, and WT_max in Eq (<ref>) are randomly generated from uniform distribution as follows:

  
  * Val_ij∼ U(1.0,30.0)    (i≠ j)
  
  * Val_i∼ U(1.0,30.0), 
  
  * WT_i∼ U(1.0,30.0), 
  
  * WT_max∼ U^integer(15,16,⋯,300),

where U^integer has the integer domain. For i=j, Val_ij=0. The penalty constant λ in Eq. (<ref>) is 100. Note that the cost function is not the QUBO formulation. There are a few ways to deal with the penalty terms for the inequality constraints. For example, we can convert the penalty terms into the QUBO formulation via duality relations <cit.>. In Ref. <cit.>, a mechanism, such as asymmetric two-body interactions and slack variables, was employed. However, we here directly evaluate the cost function value from Eq. (<ref>) for simplicity because this paper aims to check the ability of the Koopman theory. Hence, the spin number is the same as that of items; there are 15 spins.

For the evaluation of S, we use the following temperature settings:

    T_r = 10 × r + 10,

and R = 10. Other settings for the time-series data analysis are the same as in Sect. <ref>. In the solution search procedure, we use the same temperature parameters with Eq. (<ref>). The number of iterations for each search is 200, and the other settings for the search procedure are the same as in Sect. <ref>. 





Figure <ref>(a) shows the temperature-dependency of the the sum of Koopman eigenvalues. The values of R' are as follows:


  * R'=4 for Cases 1, 2, 3, 4, 7, and 9,

  * R'=5 for Case 5,

  * R'=6 for Cases 6 and 8.

Figure <ref>(b) show the relationship between the number of finding times of the optimal solution in 100 trials and the feature S; as in Sect. <ref>, there is a clear negative correlation even in the 0-1 knapsack problem.





 §.§ Data analysis with SVD


The SVD is well-known for capturing the characteristics of data matrices. Therefore, we here analyze the time-series data with the SVD. In the EDMD method, we investigated the temperature dependency of the sum of eigenvalues. By contrast, we here evaluate the temperature dependency of the sum of singular values. Instead of the quantity S in Eq. <ref>, we evaluate the following quantity:

    S^svd = ∑^R'_r=1(1/2(w^svd_r+1 + w^svd_r)(T_r+1-T_r)),

where w^svd_r is the sum of singular values at the r-th temperature.








Figure <ref>(a) shows the temperature-dependency of the the sum of singular values. The values of R' are 15 for all cases. Figure <ref>(b) show the relationship between the number of finding times of the optimal solution in 100 trials and the feature S^svd. Figures <ref>(a) and (b) shows the corresponding results for the 0-1 quadratic knapsack problems; here, R' in Eq. (<ref>) are as follows:


  * R'=5 for Case 3, 4, 6, 7, 8, and 9,

  * R'=6 for Case 2 and 5,

  * R'=8 for Case 1.

The observation from Figs. <ref>(a) and (b) are the same as the random problem cases.

From the comparison of the results of the EDMD and SVD analyses, we obtain the following observations:


  * The EDMD and SVD analyses show negative correlations between the proposed features and difficulty in problem-solving.

  * The sums of Koopman eigenvalues take small values at low temperatures in all cases, but those of singular values take various values at low temperatures.

  * While Figs. <ref>(a) and <ref>(a) show slow temperature-dependency, the EDMD results show a steep temperature dependency, as shown in Figs. <ref>(a) and <ref>(a). Hence, a narrower temperature range is sufficient for the EDMD analysis than the SVD, contributing to lower computational costs.


We consider that the steeper temperature dependency in the Koopman analysis stems from the fact that EDMD represents the dynamics itself; the SVD focuses on the entire time series data, while the EDMD deals with snapshot pairs. As the temperature increases, the dynamics should change, and the EDMD could immediately reflect the details of the dynamics because of the analysis of the snapshot pairs. Therefore, the change in dynamics with a small temperature change would be apparent in the EDMD analysis.

These numerical experiments confirm that the data analysis flow based on the Koopman analysis works well in discussing the difficulty in problem-solving.





§ CONCLUSION



In this paper, we clarified the relationship between the Koopman analysis and the difficulty in problem-solving for combinatorial optimization problems. Since only a single time-series data at a temperature is not enough to discuss the difficulty, we generated multiple time-series data at various temperatures by using the replica exchange Monte Carlo method; as denoted above, a type of specialized hardware has employed the replica exchange Monte Carlo method for the parallel computation. In the data analysis flow, we proposed the quantity that aggregates the information of the multiple time-series data. The numerical experiments confirmed the validity of the proposals; there is a negative correlation between the proposed quantity and the ability of the solution search.

This paper is the first work to use Koopman analysis to discuss the difficulty in problem-solving for combinatorial optimization problems. There are some remaining tasks left for the future. For example, one should check the algorithm dependency of the analytical results; although the Metropolis rule employed in this paper is popular in optimization, the time-series data depends on the update algorithm. It would also be interesting to evaluate the algorithm with a mechanism to get out of the local solution, for example, as proposed in Ref. <cit.>. One of the other remaining tasks is the reduction of the computational costs in the EDMD; for example, the online EDMD <cit.> could be effective for sequential data acquisition in the replica exchange Monte Carlo method. There are many recent studies for specialized hardware for combinatorial optimization problems, and customizing the algorithm for each problem would be beneficial. We believe that this paper will lead to future studies to quantify the difficulty in problem-solving.



This work was supported by JST FOREST Program (Grant Number JPMJFR216K, Japan).




10



Kirkpatrick1983
S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, Science 220, 671 (1983).




Bunyk2014
P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, and J. Whittaker, 
IEEE Trans. Appl. Supercond. 24, 1700110 (2014).


d-wave
C. McGeoch and P. Farré, D-Wave Tech. Report Series 14-1049A-A (2020).


d-wave-company
D-Wave Systems (accessed July 31, 2024). [Online] 

<https://www.dwavesys.com>.

Kadowaki1998
T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998).


Farhi2001
E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).




Takemoto2019
T. Takemoto, M. Hayashi, C. Yoshimura, and M. Yamaoka, IEEE International Solid-State Circuits Conference (ISSCC), 2019, p. 52.


Takemoto2021
T. Takemoto, K. Yamamoto, C. Yoshimura, M. Hayashi, M. Tada, H. Saito, M. Mashimo, and M. Yamaoka, IEEE Int. Solid-State Circ. Conf. 64, 64 (2021).





Goto2019a
H. Goto, K. Tatsumura, and A. Dixon, Sci. Adv. 5, eaav2372 (2019).


Goto2019b
H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura, Sci. Adv. 7, eabe7953 (2021).


Tatsumura2019
K. Tatsumura, A. R. Dixon, and H. Goto, 29th International Conference on Field Programmable Logic and Applications (FPL), 2019, p. 59.




Aramon2019
M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. G. Katzgraber, Front. Phys. 7, 48 (2019).


Matsubara2020 
S. Matsubara, M. Takatsu, T. Miyazawa, T. Shibasaki, Y. Watanabe, K. Takemoto, and H. Tamura, 25th Asia and South Pacific Design Automation Conference (ASP-DAC), 2020, p. 667.





Swendsen1986
R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 57, 2607 (1986).


Hukushima1996
K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).





Gallo1980
G. Gallo, P. L. Hammer, and B. Simeone, Combinatorial Optimization 12, 132 (1980).



Geman1984
S. Geman and D. Geman, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6, 721 (1984).
 

Suzuki2005
S. Suzuki and M. Okada, J. Phys. Soc. Jpn. 74, 1649 (2005).


Koopman1931
B. O. Koopman, Proc. Nat. Acad. Sci. 17, 315 (1931).



Rowley2009
C. Rowley, I. Mezić, P. Bagheri, S. Schlatter, and D. Henningson, J. Fluid Mech. 641, 115 (2009).



Williams2015
M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, J. Nonlinear Sci. 25, 1307 (2015).




Mauroy2020
A. Mauroy and J. M. Gonçalves, IEEE Tran. Auto. Control 65, 2550 (2020).


Korda2018
M. Korda and I. Mezić, Automatica 93, 149 (2018).



Lucas2014
A. Lucas, Front. Phys. 2, 2 (2014).





Brunton2022
S. L. Brunton, M. Budišić, E. Kaiser, and J. N. Kutz, SIAM Rev. 64, 229 (2022).


Yin2023
F. Yin, H. Tamura, Y. Furue, M. Konoshima, K. Kanda, and Y. Watanabe, J. Phys. Soc. Jpn. 92, 034002 (2023).



Sato2019
G. Sato, M. Konoshima, T. Ohwa, H. Tamura, and J. Ohkubo, Phys. Rev. E 99, 042106 (2019).



Sato2024
Y. Sato, M. Konoshima, H. Tamura, and J. Ohkubo, J. Phys. Soc. Jpn. 93, 044802 (2024).


Zhang2019
H. Zhang, C. W. Rowley, E. A. Deem, and L. N. Cattafesta, SIAM J. Appl. Dyn. Syst. 18, 1586 (2019).


























 























































