Open access peer-reviewed chapter

Optimization by Use of Nature in Physics Beyond Classical Simulated Annealing

By Masayuki Ohzeki

Submitted: November 24th 2011Reviewed: June 12th 2012Published: August 29th 2012

DOI: 10.5772/50636

Downloaded: 852

1. Introduction

We prefer to find the most appropriate choice in daily life for convenience and efficiency. When we go to a destination, we often use a searching program to find the fastest way, the minimum-length path, or most-reasonable one in cost. In such a searching problem, we mathematically design our benefit as a multivariable function (cost function) depending on many candidates and intend to maximize it. Such a mathematical issue is called the optimization problem. Simulated annealing (SA) is one of the generic solvers for the optimization problem [14]. We design the lowest-energy state in a physical system, which corresponds to the minimizer/maximizer of the cost function. The cost function to describe the instantaneous energy of the system is called as the HamiltonianH0(σ1,σ2,,σN), where σiis the degrees of freedom in the system and Nis the number of components related with the problem size. The typical instance of the Hamiltonian is a form of the spin glass, which is the disordered magnetic material, since most of the optimization problems with discrete variables can be rewritten in terms of such a physical system,

H0(σ1,σ2,,σN)=-ijJijσiσj,E1

where σiindicates the direction of the spin located at the site iin the magnetic material asσi=±1. The summation is taken over all the connected bonds (ij)through the interactionJij. The configuration of Jijdepends on the details of the optimization problem.

Then we introduce an artificial design of stochastic dynamics governed by the master equation.

ddtP(σ;t)=σ'M(σ|σ';t)P(σ';t),E2

where P(σ;t)is the probability with a specific configuration of σisimply denoted as σat timet. The transition matrix is written asM(σ'|σ;t), which obeys the conservation of probability σM(σ|σ';t)=1and the detailed balance condition

M(σ|σ';t)Peq(σ';t)=M(σ'|σ;t)Peq(σ;t).E3

Here we denote the instantaneous equilibrium distribution (Gibbs-Boltzmann distribution) as

Peq(σ;t)=exp(-β(t)E(σ;t))/Zt,E4

where the instantaneous energy E(σ;t)is the value of the Hamiltonian H0and Ztdenotes a normalization factor termed as the partition function. In order to satisfy these conditions, we often use the transition matrix with Metropolis rule as

M(σ|σ';t)=min(1,exp(-βΔE(σ|σ';t))),E5

where

ΔE(σ|σ';t)=E(σ;t)-E(σ';t),E6

or heat-bath rule as

M(σ|σ';t)=δ1(σ,σ')exp-β2ΔE(σ|σ';t)2coshβ2ΔE(σ|σ';t),E7

where

δ1(σ,σ')=δ(2,i=1N(1-σiσ'i)).E8

The master equation simulates behavior of relaxation toward a specific distribution associated with the energy of the system. If we evolve the system for a long time with a virtual parameter β(t)being a constantβ, which is the inverse temperature in context of physics, the probability distribution converges to the equilibrium distribution. To generate lower energy state, let us setβ1. Unfortunately, the spin-glass system in the low-temperature often exhibits the extremely long time for equilibration of the system. The most relevant reason is on the complicated structure of the energy of the spin-glass system as schematically depicted in Fig. 1.

Figure 1.

Energy structure of a spin-glass system. We map the possible 2N configurations to the one-dimensional horizontal axis for simplicity. The vertical axis represents the value of the energy for each configuration.

There are barriers between the valleys to avoid hopping from state to state in the low-temperature, where the energy effect is dominant. Therefore it is difficult to reach the equilibrium distribution by a direct simulation with a constant temperature. Instead, by tuning a virtual parameter β(t)from zero to a large number, we perform stochastic searching while keeping to trace the instantaneous equilibrium distribution. The mathematical guarantee to converge to the instantaneous equilibrium state by gradually changing the inverse temperature has been proved by Geman and Geman based on the classical inhomogenious (time-dependent) Markov chain representing nonequilibrium processes [6]. The convergence theorem states that we reach the equilibrium distribution, if we obey the following schedule or slower rate to change the inverse temperature as

β(t)=1pNlogαt+1,E9

where αis exponentially small in Nand pis a constant independence ofN. An intuitive way to understand the performance of SA is as follows. The inverse temperature controls the range of searching, roughly speaking. The instantaneous state keeps hopping from state to state in a relatively high-temperature region as depicted in Fig. 2.

Figure 2.

Behavior of the instantaneous state in SA. The left panel shows the stochastic searching in SA in the high-temperature regionβ≈0, while the right one describes that in the low-temperature one. The dashed curves express the structure of the energy as in Fig. 1, while the thick ones denote the probability for realization of each configuration schematically.

By gradually decrease of the inverse temperature, we narrow the range of searching. The lower energy state means its realization with a higher possibility following the instantaneous equilibrium distribution as in Fig. 2. Demand of a sufficiently slow control of the inverse temperature implies that we need enough time to find the states with relatively lower energies by the stochastic searching before the barrier avoids globally searching for the lower energy state.

Basically, SA is based on the behavior closely to the instantaneous equilibrium state. Therefore we need to perform the change of the inverse temperature with a sufficiently slow control. In order to improve the performance, in particular to shorten the necessary time, we need to consider the protocol away from the equilibrium state, that is nonequilibrium process.

In this chapter, we show a novel direction to solve efficiently the optimization problem by use of the nature in nonequilibrium physics. In statistical physics, the interest of researchers in nonequilibrium dynamical behavior has increased. Among several remarkable developments, the Jarzynski equality (JE), which is known as a generalization of the second law of thermodynamics, might be possible to change the paradigm in optimization problem by use of the physical nature. The Jarzynski equality relates an average over all realizations during a predetermined nonequilibrium process with an expectation in an equilibrium state. As seen later, the mathematical structure of JE does not depend on the schedule and the rate of changing the external parameter. It means that, if we implement JE to solve the optimization problem, we do not need to demand slow control of the driver of the system. The challenge of the implementation of JE have been performed in several researchers. Although not yet have been studied the performance in the actual application to the optimization problem, we show the possibility of the novel method from several analyses.

2. Population annealing

We introduce a couple of theories in nonequilibrium statistical physics in short, before we show the actual application. They provide the supplement to make the protocol of SA faster. The Jarzynski equality is the most important key.

2.1. Jarzynski equality

Among several recent developments in nonequilibrium statistical mechanics, we take JE as an attempt to improve the performance of SA. The Jarzynski equality relates quantities at two different thermal equilibrium states with those of nonequilibrium processes from t=0to t=Tconnecting these two states as [10,11]

e-βW0T=ZTZ0,E10

where the partition functions appearing in the ratio on the right-hand side are for the initial (t=0) and final Hamiltonians (t=T). The quantity on the right-hand side can be represented by the exponentiated difference of the free energy exp(-βΔF)between the initial and final conditions. The brackets on the left-hand side express the nonequilibrium average over all the instantaneous realizations of the degrees of freedom, for instance spin configurations, during the nonequilibrium process with the following path probability defined as

P0T({σt})=k=0n-1eδtM(σk+1|σk;tk)Peq(σ0;t0).E11

It implies that the observations of the nonequilibrium behavior can estimate the equilibrium quantity represented by the partition functions, that is the free energy. This equality is regarded as a generalization of the well-known inequality, the second law of thermodynamicsW0TΔF, which can be reproduced by the Jensen inequality. One of the important features is that JE holds independently of the pre-determined schedule of the nonequilibrium process.

In order to consider the improvement of SA, let us apply the nonequilibrium process with change of the temperature. We then have to employ the pseudo work instead of the ordinary performed work due to the energy difference as

Y(σ;tk)=βk+1-βkE(σ),E12

where we use discrete time expressions as t0=0and tn=Tfor simplicity and we assume that the instantaneous energy does not depend on the time asE(σ). The Jarzynski equality holds formally in the case with change of temperature,

e-Y0T=ZTZ0.E13

We show a simple proof of JE for the particular dynamics in SA. Let us consider a nonequilibrium process in a finite-time schedule governed by the master equation. The left-hand side of JE is written as

e-Y0T={σk}k=0n-1e-Y(σk+1;tk)eδtM(σk+1|σk;tk)Peq(σ0;t0).E14

where we use the formal solution of the master equation by the exponentiated transition matrix. We take the first product of the above equation as,

σ0e-Y(σ1;t0)eδtM0(σ1|σ0;t0)Peq(σ0;t0)E15
    =Peq(σ1;t1)Zt1Z0.E16

Repetition of the above manipulation in Eq. (14) yields the quantity in the right-hand side of JE as,

σnPeq(σn;tn)k=0n-1Ztk+1Ztk=ZTZ0.E17

2.2. Fluctuation theorem

The Jarzynski equality is a consequence from the fluctuation theorem [2, 3, 4], which relates the probability P0T({σt})with that of the inverse process PT0({σt})as

P0T({σt})PT0({σt})e-Y=ZTZ0.E18

This leads to the more generic result, for an observable O({σt})depending on the instantaneous spin configurations,

O({σt})e-Y0T=Or({σt})T0ZTZ0,E19

where Ordenotes the observable which depends on the backward processT0. The brackets with subscript 0Texpress the average with the weight P0T({σt})over possible realizations{σt}. ForO=Or=1, Eq. (18) reduces to JE.

If we choose an observable depending only on the final state, which is denoted asOT, instead ofO({σt}), Orreads an observable at the initial state in the backward process. Then OrT0equals to the ordinary thermal average at the initial equilibrium state with βTrepresented byβT, and Eq. (18) leads to

OTe-Y0T=OβTZTZ0.E20

By looking over the above calculations, we can understand the roll of the exponentiated pseudo work. The resultant distribution after SA is given byP0T({σt}). The biased distribution with the exponentiated pseudo work, P0T({σt})exp(-Y)always gives the equilibrium one. In this sense, the exponentiated pseudo work plays a roll to fill the gap between the equilibrium distribution and the resultant one after SA. Therefore, if we skillfully use the exponentiated pseudo work to keep the instantaneous distribution close to the equilibrium one, we can invent the improved version of SA.

2.3. Population annealing

We introduce an improvement of SA by use of the property of JE. Let us consider to implement JE in numerical simulation. We parallelize the instantaneous spin configurations {σt}as{σt}i=1,,C. Each of the configuration independently will be evolved by the master equation. We regard the exponentiated pseudo work on the left-hand side of JE, Eq. (14), as the weight for each realization of the configuration. By computing the pseudo work for each realization and multiplying the weight given by the exponentiated pseudo work, we simultaneously perform the stochastic dynamics governed by the master equation. At the last stage of repetition of the above procedure, we obtain the ratio of the partition function as in the right-hand side of JE. In order to calculate its value, we estimate the empirical average as, after parallel computing of the master equation,

1Ci=1Cexp-k=1nY(tk;σ).E21

While estimating the ratio of the partition functions by JE, implementation of Eq. (19) gives the thermal average of the observable through their ratio. This is the typical implementation of JE in a numerical simulation, which is called as population annealing (PA) [7, 9, 17] as depicted in Fig. 3,

Figure 3.

Schematic picture of the process of PA byC=4. The size of the circles denotes the weight given by the multiplication of the exponentiated pseudo work during PA.

We remark that, as proposed in the literatures [7, 9], we have to employ a skillful technique, resampling, to efficiently generate the relevant copies to estimate the nonequilibrium average and maintain the stability of the method. The population annealing with resampling method indeed shows outstanding performance comparable to a standard technique to equilibrate the spin-glass system known as the exchange Monte Carlo method [8]. If we successfully generate the equilibrium distribution in the low-temperature region, we efficiently find the lowest energy state, which corresponds to the optimal solution in context of the optimization problem. Therefore PA is also relevant for the improvement of SA as a solver for the optimization problem. The advantage of PA is cutting the computational time compared to SA, since PA follows the property of JE. It means that we find the optimal solution by use of PA faster than SA.

Below, we propose an ambitious use of PA to evaluate the equilibrium property in the low-temperature in spin glasses by use of the special symmetry [21].

2.4. Spin glass

At first we briefly review a useful analysis in several spin-glass systems, which provides a powerful technique to discuss the possibility of PA. Let us consider a simple model of spin glasses, the ±JIsing model, on an arbitrary lattice. The Hamiltonian is the same form as Eq. (1) as

H=-ijJτijσiσj,E22

where we extract the sign of the interaction τijfollowing the distribution function as

P(τij)=pδ(τij-1)+(1-p)δ(τij+1).E23

The partition function, which is the most important quantity through the free energy, is defined as

Z(K;{τij})={σi}ijexp(Kτijσiσj).E24

The free energy is then given by

-βF(K;{τij})=logZ(K;{τij}),E25

where the product βJis rewritten asK. Both of the above quantities depend on the specific configuration{τij}. In order to evaluate the physical property of spin glasses in equilibrium, we strive the difficult task to deal with the free energy depending on the non-uniform interactions. Instead of the direct manipulation, the averaged quantity over all the possible configurations of {τij}may be considered based on the self-averaging property as, in the large-limitN,

1NF(K;{τij})1NF(K;{τij}),E26

where the square bracket denotes the average over all the combinations of {τij}(configurational average). The self-averaging property is valid for other observables, which can be obtained from the free energy per site like the internal energy.

2.5. Gauge transformation

Here let us define a local transformation by the simultaneous change of the interactions and spin variables as, by the binary variables εi=±1[18, 19]

τijεiεjτijE27
σiεiσi.E28

This is called as the gauge transformation. Notice that the gauge transformation does not alter the value of the physical quantity given by the double average over τijand σisince it changes only the order of the summations. The Hamiltonian can not change its form after the gauge transformation since the right-hand side is evaluated as

-ijJτijεiεjεiσiεjσj=H,E29

As this case, if the physical quantity is invariant under the gauge transformation (gauge invariant), we can evaluate its exact value even for finite-dimensional spin glasses. The key point of the analyses by the gauge transformation is on the form of the distribution function. Before performing the gauge transformation, the distribution function can take the following form as

P(τij)=eKpτij2coshKp,E30

whereexp(-2Kp)=(1-p)/p. The gauge transformation changes this into

P(τij)=eKpτijεiεj2coshKp.E31

Let us evaluate the internal energy by aid of the gauge transformation here. The thermal average of the Hamiltonian is given by

HK={σi}1Z(K;{τij})HijexpKτijσiσjE32
=-JddKlogZ(K;{τij}).E33

We can use the self-averaging property here and thus take the configurational average as

HKKp={τij}ijexp(Kpτij)2coshKp×HK.E34

where []Kpdenotes the configurational average by the distribution function (29) withKp. Then we perform the gauge transformation, which does not change the value of the internal energy due to gauge invariance,

HKKp={τij}ijexp(Kpτijεiεj)2coshKp×HK.E35

Therefore we here take the summation over all the possible configurations of {σi}and divide it by 2N(the number of configurations) as

HKKp=12N{εi}{τij}ijexp(Kpτijεiεj)2coshKp×HK.E36

We take the summation over {εi}in advance of that over {τij}and then find the partition function with Kpinstead ofK.

HKKp=12N{τij}Z(Kp;{τij})(2coshKp)NB×HK.E37

where NBis the number of bonds. Going back to the definition (31), both of the partition functions on the denominator and numerator can be cancelled when Kp=Kas

HKK=-J2N(2coshKp)NB{σi}{τij}ddKexpKτijσiσjE38
=-NBJtanhK.E39

Similarly, we can evaluate the rigorous upper bound on the specific heat. The condition Kp=Kconfines the special subspace in which we can perform the exact analysis for spin glasses. This subspace is called as the Nishimori line (NL) [18, 19].

2.6. Jarzynski equality for spin glasses

By use of the gauge transformation as above introduced briefly, let us consider the application of the relations (18) and (19) to spin glasses, namely PA for such a complicated system in a tricky way. We analyze JE for the spin-glass model for several interesting quantities below.

2.6.1. Gauge-invariant quantities like internal energy

We apply Eq. (19) to a gauge-invariant quantity G({τij})for the purpose of evaluation of the equilibrium quantity in spin glasses. The configurational average for Eq. (19) yields

GT({τij})e-YK0KTKp=G({τij})KTZ(KT;{τij})Z(K0;{τij})Kp.E40

The quantity on the left-hand side is the configurational and nonequilibrium averages of the observable G({τij})at final stage of PA, that is after the protocol K0KTwith the factore-βW. On the other hand, G({τij})KTon the right-hand side expresses the configurational and thermal averages of the equilibrium state for the final Hamiltonian.

The gauge transformation σiσiεi,τijτijεiεj(i,j)leads us to

G({τij})KTZ(KT;{τij})Z(K0;{τij})Kp={τij}G({τij})KTijeKpτijεiεj(2coshKp)NBZ(KT;{τij})Z(K0;{τij}).E41

All the quantities in this equation are invariant under the gauge transformation. The summation over {εi}and division by 2Ngives

G({τij})KTZ(KT;{τij})Z(K0;{τij})Kp={τij}G({τij})KTZ(Kp;{τij})2N(2coshKp)NBZ(KT;{τij})Z(K0;{τij}).E42

On the other hand, let us evaluate the quantityG({τij})KTKp. Similarly to the above calculation, the following identity can be obtained by the gauge transformation,

G({τij})KTKp={τij}G({τij})KTZ(Kp;{τij})2N(2coshKp)NB.E43

By Setting Kp=K0in Eq. (40) and Kp=KTin the above equation, we reach the following nonequilibrium relation,

GT({τij})e-YK0KTK0=G({τij})KTKTcoshKTcoshK0NB.E44

If we set GT({τij})=1in the resultant equation, the Jarzynski equality for spin glass is obtained,

e-YK0KTK0=coshKTcoshK0NB.E45

Equation (43) leads to the lower bound on the pseudo work, using Jensen’s inequality for the average ofe-Y,

YK0KTK0-NBlogcoshKTcoshK0.E46

By substituting GT({τij})=Hinto Eq. (42), we obtain

He-βWK0KTK0=HKTKTcoshKTcoshK0NB.E47

This equation shows that the internal energy after the cooling as in SA or heating process starting from a temperature on NL, which is in the present caseK0=Kp, is proportional to the internal energy in equilibrium on NL corresponding to the final temperature KT=Kpas in Fig. 4. We here assume to perform PA and take an average over all results after many repetitions. The nonequilibrium process starts from NL (1/K0,1/K0)and ends at the point away from NL(1/K0,1/KT). Notice that the ordinary procedure in PA gives the estimation of the equilibrium quantity at the last condition as(1/K0,1/KT). However the corresponding internal energy is at the different point but on NL as(1/KT,1/KT). It means that we can obtain the equilibrium quantities in the different amount of the randomness from the initial condition through the configurational average of the results from PA.

Figure 4.

The nonequilibrium process as in SA/PA on the left-hand side as in Eq. (45) drawn as an arrow on the phase diagram. The corresponding equilibrium state as on the right-hand side of Eq. (45) is at the point on NL. The solid curves are the phase boundaries and the dashed line of 45o represents NL.

2.6.2. Gauge-non-invariant quantities

In statistical physics, it is important to detect the order of the instantaneous spin configuration in the system. For instance, as in Fig. 4, there are several phases, ferromagnetic, paramagnetic and spin-glass ones, involved in the spin-glass model. They have the characteristic quantities to distinguish themselves, termed as the order parameter. The order parameter to identify the phase boundary between the ferromagnetic and paramagnetic phases is the magnetization defined as

m=1Ni=1NσiE48

Therefore it is important to observe the behavior of the first momentum of spin variable in equilibrium. For this purpose, we choose σi(T)for Oin Eq. (19) and consider the possibility of the application of PA. After the configurational average, we obtain

σi(T)e-YK0KTKp=σiKTZ(KT;{τij})Z(K0;{τij})Kp.E49

Gauge transformation for the right-hand side in this equation yields

σiKTZ(KT;{τij})Z(K0;{τij})Kp={τij}Z(KT;{τij})Z(K0;{τij})σiKTεiijeKpτijεiεj2coshKp.E50

We again sum both sides of this equation over all the possible configurations of {εi}and divide the obtained quantity by 2Nto find

σiKTZ(KT;{τij})Z(K0;{τij})Kp={τij}Z(Kp;{τij})2N(2coshKp)NBσiKTεiKpZ(KT;{τij})Z(K0;{τij}).E51

The following relation can also be obtained in a similar manipulation,

σiK0Kp={τij}Z(Kp;{τij})2N(2coshKp)NBσiK0εiKp.E52

By setting Kp=K0in Eq. (49) and Kp=KTin Eq. (50), we reach a relation

σiKTZ(KT;{τij})Z(K0;{τij})K0=σiK0KTcoshKTcoshK0NB.E53

As a result, we obtain a nonequilibrium relation,

σi(T)e-YK0KTK0=σiK0KTcoshKTcoshK0NB.E54

The same method yields another relation for the correlation functions to similarly measure the magnitude of order in the system

σ0(T)σr(T)e-YK0KTK0=σ0σrK0KTcoshKTcoshK0NB.E55

The obtained relations (52) and (53) relate the equilibrium physical quantities away from NL (the right-hand sides) with other quantities measured during the nonequilibrium process from a point on NL to another point away from NL (the left-hand sides) as depicted in Fig. 5. The spin-glass system in the low-temperature region exhibits the extremely slow relaxation toward equilibrium. This feature hampers to observe the equilibrium behavior of spin glasses. However our results imply that the configurational average of PA would overcome the difficulty. One may attempts the heating process from NL in order to evaluate the low-temperature property through Eqs. (52) and (53) as depicted in Fig. 5. The Jarzynski equality holds irrespectively of the schedule to control the external field. It means that we can investigate the low-temperature behavior for spin glasses without suffering from critical slowing down.

Figure 5.

Equations (52) and (53) relate equilibrium quantities at the point (1/K0,1/KT) (the lower-left dot) with other physical quantities estimated during the nonequilibrium process shown in the arrow. The lower-left dot is in the spin glass phase whereas the corresponding arrow is in the ferromagnetic phase.

Unfortunately, however, the exponentiated pseudo work does not hold the self-averaging property. It means that the sample-to-sample fluctuation between different configurations of {τij}remains to be relevant even in a large-Nsystem. Therefore, if we estimate the empirical average of realizations of {τij}following the obtained equalities, we do not correctly reproduce the quantity of the right-hand side. We describe the test results for estimation of the ratio of the partition functions as in Eq. (13) by PA for 100realizations of {τij}in Fig. 6. We perform PA for the ±JIsing model on the square lattice with a linear sizeL=6, and Monte-Carlo stepR=1000. The number of copies isC=100.

Figure 6.

Superimposed plots of Differences between the exact value and the estimation given by PA. The horizontal axis denotes the instantaneous temperature during PA. The vertical axis represents the difference from the exact values given by the transfer matrix method in advance.

The population annealing can correctly estimate the ratio of the partition functions for each realization, but their simple average does not coincide with the quantity on the right-hand side of Eq. (43) as in Fig. 7. Both of the results are away from the exact solutions due to the sample-to-sample fluctuation and show nontrivial behavior depending on the linear size. These facts imply lack of self-averaging property. Therefore, if we exploit all the above results given by the configurational average of the exponentiated pseudo work, we have to overcome this violation due to lack of the self-averaging property. This is one of the remaining problem associated with this procedure with PA for spin glasses.

Figure 7.

Difference between the exact values and the empirical averages over 100 realizations. The cross points give the deviation of the empirical averages of the exact values of ratios of the partition functions from the exact value on the right-hand side of Eq. (43), while the tilted-cross ones represent that of the empirical average of the estimation given by PA.

3. Quantum annealing

Observant readers may begin to recognize the possibility to use physical nature to drive the system in searching the lowest energy (ground state) instead of thermal fluctuation controlled by the inverse temperature. We show another strategy to find the ground state recently studied in a field of physics, quantum annealing (QA) [13].

3.1. Quantum adiabatic computation

In quantum-mechanical system, we can use a parallel way to drive all the candidates of the desired solution in optimization problem by use of superposition. Quantum annealing uses quantum fluctuation between superposed states to search for the ground state. One of the successful strategies is to use the adiabatic evolution known as quantum adiabatic computation (QAC) [5]. In QAC, as the procedure of SA, we control to gradually decrease the strength of quantum fluctuations to drive the system. Similarly to SA, the convergence into the optimal solution of QAC (the ground state) is also guaranteed by a mathematical proof [15].

In QAC, we introduce a non-commutative operator to drive the system by quantum nature in addition to the original HamiltonianH0, which is designed to represent the optimization problem to be solved, as

H(t)=f(t)H0+1-f(t)H1,E56

where f(t)is assumed to be a monotonically increasing function satisfying f(0)=0andf(T)=1. For instance, f(t)=t/T, where Tdenotes the computation time for QAC. In order to exemplify the explicit instance, we again assume that H0is the spin glass Hamiltonian as

H0=-ijJijσizσjz,E57

where σizis the zcomponent of the Pauli operators defined as

σx=0110,    σy=0i-i0,    σz=100-1.E58

We take the computational basis of the eigenstates of the z-component of the Pauli matrix (Ising variables) to represent the instantaneous state as|Ψ(t)=|σ1z,σ2z,,σNz. The transverse-field operator is often used as quantum fluctuations for implementing QAC for the spin-glass model

H1=-Γ0i=1Nσix.E59

where Γ0is the strength of the transverse field. The whole Hamiltonian of QAC (although widely used also for QA) thus becomes

H(t)=f(t)ijJijσizσjz+1-f(t)Γ0i=1Nσix.E60

The quantum adiabatic computation starts from a trivial ground state ofH1. In the present case, the ground state of the transverse-field operator H1is simply written by a uniform linear combination as{σ}|σ1z,σ2z,,σNz/2N.

The adiabatic theorem guarantees that the instantaneous state at timet, |Ψ(t), is very close to the instantaneous ground state for a sufficiently large T(implying slow control) as|0(t), if the instantaneous ground state |0(t)is non-degenerate. The condition for|0(t), 0(t)|Ψ(t)1-δ2(δ1)to hold can be written as [15]

max1(t)|dH(t)dt|0(t)minΔ2(t)=δ,E61

where |1(t)is the instantaneous first excited state, and Δ(t)is the energy gap between the ground state and first excited one. The maximum and minimum should be evaluated between 0andT. In the present case, sincedH(t)/dt1/T, the above adiabatic condition is reduced into

T1δminΔ2(t).E62

It means that if we desire to solve the optimization problems by use of QAC, which one of the specialized version of QA, we take the computational time proportional to the inverse square of the energy gap. If the problem involved with the exponential closure of the energy gap for increasing ofN, QAC must take extremely long time to obtain the ground state with a high probability [12, 25]. Interestingly, we can reproduce the convergence theorem of SA (9) from the above adiabatic condition with recourse to a mathematical mapping of the procedure of SA into the quantum dynamics [15, 23]. It implies that the nature of QAC can be understood through that of SA by the mathematical mapping technique.

Below, we would provide a new paradigm to solve faster than the ordinary scheme of SA. A fast sweep of the system yields nonequilibrium behavior. Although we have not yet understood deeply the nonequilibrium phenomena, there are a few well-established theories which rises to applications to the optimization problem. One possibility is PA for the quantum system by use of JE and its alternatives [20]. Here we again employ JE to give another scheme of QA while considering the nonequilibrium behavior.

3.2. Jarzynski equality for isolated quantum system

In order to consider the nonequilibrium behavior away from the adiabatic dynamics of QAC, we shortly review JE for an isolated quantum system [1, 24].

To directly use JE in the protocol to find the ground state of the spin-glass Hamiltonian H0as in Eq. (55), we prepare the dynamical quantum system following the time-dependent Hamiltonian (58). In addition, we pick up a state from the canonical ensemble for H(0)=H1=-Γ0iσixat the initial stage of the procedure and then let it evolve following the time-dependent Schrödinger equation. We measure the performed work in the isolated quantum system asW=Em(T)-En(0), which is simply given by the difference between the outputs of projective measurements of the initial and final energies. Here mand ndenote the indices of the instantaneous eigenstates as H(T)|m(T)=Em(T)|m(T)andH(0)|n(0)=En(0)|n(0), respectively. The time-evolution operator is given by the following unitary operator as

UT=Texpi0TdtH(t),E63

where Tdenotes the time ordered product. Therefore the transition probability between the initial and final stages is given as

Pm,n(0T)=|Ψm(T)|UT|Ψn(0)|2.E64

The path probability for the nonequilibrium process starting from the equilibrium ensemble is then evaluated as

Pm,n(0T)exp(-βEn(0))Z0(β;{Jij}),E65

where we express the instantaneous partition function for the instantaneous Hamiltonian at each time tasZt(β;{Jij}).

By directly evaluating the nonequilibrium average of the exponentiated work but for the isolated system, we reach JE applicable to non-adiabatic version of QA. We define the nonequilibrium average of the exponentiated work as

e-βWQA=m,ne-βWPm,n(0T)exp(-βEn(0))Z0(β;{Jij}),E66

which becomes the left-hand side of JE. The quantity defined here is evaluated as

e-βWQA=m,ne-βEm(T)Z0(β;{Jij})Pm,n(0T)E67
=me-βEm(T)Z0(β;{Jij})E68
=ZT(β;{Jij})Z0(β;{Jij}),E69

where we used the fact that the performed work Wwas a classical number and

nPm,n(0T)=nΨm(T)|UT|Ψn(0)Ψn(0)|UT|Ψi(T)E70
=mΨm(T)|UTUT|Ψm(T)=1.E71

If we measure the physical observable O^Tat the last of the nonequilibrium process, we obtain another equation as, similarly to the classical version,

O^Te-βWQA=O^βZT(β;{Jij})Z0(β;{Jij}),E72

where the subscript on the square brackets in the right-hand side denotes the thermal average in the last equilibrium state with the inverse temperatureβ. The ratio of Eqs. (64) and (66) gives

O^Te-βWQAe-βWQA=O^β.E73

The resultant equation suggests that we can estimate the thermal average under the Hamiltonian H0on the right-hand side through the left-hand side of JE. This fact may be useful in the evaluation of equilibrium average, since the left-hand side is evaluated without slow adiabatic processes. Differently from QAC, we must sweep the quantum system repeatedly to correctly estimate the nonequilibrium average as in JE, but in short time. We propose such a procedure as the alternative of QAC, the non-adiabatic quantum annealing (NQA) based on the property of JE as above established. First we discuss the possibility as a solver of the optimization problem below.

4. Non-adiabatic quantum annealing

In order to investigate the property of the ground state, we tune the inverse temperature into a very large valueβ1. The nonequilibrium average on the left-hand side of JE involves a non-extensive quantity, the exponentiated work, whose value fluctuates significantly from process to process. Therefore the average on the left-hand side must be calculated by many trials of annealing processes. Thus, rare events with large values of the exponentiated work (i.e.β|W|Γ0) would contribute to the average significantly, and we have to repeat the annealing process very many, typically exponentially, times in order to reach the correct value of the average. This property would be a bottleneck of the simple implementation of NQA, instead of long time involved by the closure of the energy gap in QAC. What about PA in the classical counter part? In order to generate the relevant contributions in PA, we use the biased distribution with the exponentiated pseudo work through resampling technique [7, 9]. Without resampling technique, we cannot efficiently reproduce the prediction given by JE. In this fact, if we implement the biased distribution in the quantum system, we would use NQA without suffering from rare events. It means that we can perform NQA in order to find the ground state in a short time with several repetitions. So far, it is the very difficult task to realize the biased distribution in the quantum system. However it is worthwhile to consider its possibility in the future.

4.1. Several analyses of non-adiabatic quantum annealing

Unfortunately, we have not reached any positive answers on the performance of NQA. Instead let us here evaluate several properties in nonequilibrium process as in NQA for the particular spin glasses. We can exactly analyze nonequilibrium behavior by combination of JE with the gauge transformation, although there are few exact results in nonequilibrium quantum dynamical system with many components [22].

Following the prescription of JE, let us consider a repetition of NQA starting from the equilibrium ensemble. The initial Hamiltonian of NQA is given only by the transverse fieldH(0)=H1. It turns out that the starting point of our analyses is the specialized JE to the case for NQA as

e-βWQA=ZT(β,{Jij})(2coshβΓ0)N.E74

We assume that the interactions {Jij}follow the distribution function for the ±JIsing model (22), which is better to be rewritten as

P(Jij)=exp(βpJij)2coshβpJ,E75

where we do not use K=βJfor transparency, andexp(-2βpJ)=(1-p)/p.

4.2. Gauge transformation for quantum spin systems

For several special spin glasses as the ±JIsing model, the gauge transformation is available for analyses on the dynamical property even under quantum fluctuations. The time-dependent Hamiltonian as in Eq. (58) is invariant under the following local transformation,

σixσix,σiyξiσiy,σizξiσiz,JijJijξiξj    (i,j),E76

where ξi(=±1)is called as a gauge variable. This transformation is designed to preserve the commutation relations between different components of Pauli matrices [16].

4.3. Relationship between two different paths of NQA

Below, we reveal several properties inherent in NQA by the gauge transformation. Let us take the configurational average of Eq. (68) over all the realizations of {Jij}for the special case with β=β1and βp=β2as

e-β1WQAβ2=ZT(β1;{Jij})2coshβ1Γ0Nβ2.E77

The right-hand side is written explicitly as

e-β1WQAβ2={Jij}exp(β2ijJij)(2coshβ2J)NBZT(β1;{Jij})2coshβ1Γ0N.E78

Let us here apply the gauge transformation as introduced above. Since the time-dependent Hamiltonian is invariant, we may sum over all possible configurations of the gauge variables {ξi}and divide the result by 2Nin order to obtain the quantity on the left-hand side,

e-β1WQAβ2={Jij}ZT(β2;{Jij})ZT(β1;{Jij})2N(2coshβ2J)NB2coshβ1Γ0N.E79

A similar quantity of the average of the exponentiated work for spin glass with the inverse temperature β2and the parameter for the quenched randomness β1gives

e-β2WQAβ1={Jij}ZT(β2;{Jij})ZT(β1;{Jij})2N(2coshβ1J)NB2coshβ2Γ0N.E80

Comparing Eqs. (73) and (74), we find the following relation between two different non-adiabatic processes,

e-β1WQAβ2=e-β2WQAβ1coshβ1Jcoshβ2JNBcoshβ2Γ0coshβ1Γ0N.E81

We describe the two different paths of NQA related by this equality in Fig. 8.

Figure 8.

Two different processes of NQA in Eq. (75). The left-hand side of Eq. (75) represents the process toward the upper-right red dot and the right-hand side ends at the lower-left dot. Three phases expressed by the same symbols as in Fig. 4 are separated by solid curves and the dotted line describes NL (βp=β).

Setting β2=0in Eq. (75), (implyingp=1/2, the symmetric distribution), we find a nontrivial relation on the performed work during NQA

e-β1WQA0=(coshβ1J)NB(coshβ1Γ0)N.E82

The symmetric distribution (β2=0on the left-hand side) makes it possible to reduce the right-hand side to the above trivial expression. It is remarkable that NQA, which involves very complex dynamics, satisfies such a simple identity irrespective of the speed of annealingT. The Jensen inequality for the above equality leads us to the lower bound for the performed work as

WQA0-Nβlog(coshβJ)zcoshβΓ0.E83

where zis the coordination number asz=NB/N. Here we generalize the inverse temperature to βfrom the specific choiceβ1. This lower bound is loose, since the direct application of the Jensen inequality to JE for NQA yields, after the configurational average with the symmetric distribution,

WQA01βD(0|β)-Nβlog(coshβJ)zcoshβΓ0.E84

where D(β|β')is the Kullback-Leibler divergence defined as

D(β|β')={Jij}P~β({Jij})logP~β'({Jij})P~β({Jij})E85

Here we defined the marginal distribution for the specific configuration {Jij}summed over all the possible gauge transformations,

P~β({Jij})=12N{ξi}ijP(Jij)=ZT(β;{Jij})2N(2coshβJ)NB.E86

Since the Kullback-Leibler divergence does not become non-negative, the work performed by the transverse field during a nonequilibrium process in the symmetric distribution (i.e. the left-hand side of Eq. (78)) does not lower below the second quantity on the right-hand side of Eq. (78), namely Eq. (77). This fact means that Eq. (77) is a loose lower bound.

4.4. Exact relations involving inverse statistics

Beyond the above results, we can perform further non-trivial analyses for the nonequilibrium process in the special conditions. Let us next take the configurational average of the inverse of JE, Eq. (68), as

1e-βWQAβp=2coshβΓ0NZT(β;{Jij})βp.E87

The application of the gauge transformation to the right-hand side yields

1e-βWQAβp={Jij}exp(βpijJijξiξj)(2coshβpJ)NB2coshβΓ0NZT(β;{Jij}).E88

Summation of the right-hand side over all the possible configurations of {ξi}and division of the result by 2Ngive

1e-βWQAβp={Jij}ZT(βp;{Jij})2N(2coshβpJ)NB2coshβΓ0NZT(β;{Jij}).E89

This equation reduces to, by settingβp=β, namely on the Nishimori line,

1e-βWQAβ=(coshβΓ0)N(coshβJ)NB.E90

Comparison of Eqs. (76) and (84) gives

e-βWQA0=1e-βWQAβ-1.E91

As shown in Fig. 9, the resultant equation leads us to a fascinating relationship of the two completely different processes through the inverse statistics. One denotes NQA toward the Nishimori line, while the other expresses for the symmetric distribution.

Figure 9.

Two different nonequilibrium processes in NQA through Eq. (85). The same symbols are depicted as in Fig. 4. The blue circle represents the determination of the process on the right-hand side of Eq. (85), whereas the red one is for the left-hand side.

We can find the exact results through the inverse statics of the inverse statics of Eq. (66). Let us further consider the case for the two-point correlationOT=σizσjz. Taking the configurational average of both sides under the conditionβp=β, we find

1σizσjze-βWQAβ=(coshβΓ0)N(coshβJ)NB1σizσjzββ.E92

The quantity on the right-hand side becomes unity by the analysis with the gauge transformation as has been shown in the literatures [18, 19]. We thus reach a simple exact relation

1σizσjze-βWQAβ=(coshβΓ0)N(coshβJ)NB,E93

which is another exact identity for processes of NQA.

We obtain several exact nontrivial relations between completely different paths in NQA as shown above by use of the gauge transformation, which is a specialized tool to analyze spin glasses. We should notice that such results are very rare for the nonequilibrium behavior in disordered quantum system. The importance of the above equalities is still not clear. We emphasize that, when we realize the quantum spin systems in experiments, the above results would be a valuable platform to confirm their precisions and conditions. The quantum annealing is originally intended to be a tool implemented in so-called quantum computer. Therefore we need theoretical studies not only on the performance in the application but also how to make good condition to implement the protocol. In this regard, our studies shed light on the future indicator to open the way to realize the generic solver in the quantum computer. We must continue the active studies in this direction. 1.04

5. Conclusion

We reviewed several recent active studies on the alternatives of SA with use of the novel substantial progress in statistical physics. The key point was to exploit the nonequilibrium behavior during performing the active control on the system. The Jarzynski equality states the possibility to estimate the equilibrium quantities by the average quantity through the nonequilibrium behavior. It means that we can invent several new strategies by use of JE away from the paradigm of the simulated annealing, which sticks to the quasi-static dynamics. The population annealing is a starting point of the studies in this direction. It is certain that population annealing find out the desired solution of the optimization problem from the property of JE faster than SA. Roughly speaking, the population annealing cuts the computational time (CPU) by use of the parallel dynamics (memory). The remaining problem is to evaluate its qualitative performance of PA for the optimization problem.

Not only the direct use of PA, we propose another type of its application in this chapter. Regarding on this type, we show skillful analyses by use of the special symmetry hidden in spin glasses to give several nontrivial exact relations. The resultant relations are useful to investigate the low-temperature region for spin glasses if we implement them by aid of PA, since we do not suffer from the critical slowing down peculiar in spin glass.

Meanwhile, if we employ a different rule to drive the system, we would be able to find the way to solve the optimization problem as SA. We reviewed QA, which was by use of quantum fluctuations as a driver. The specialized version of QA, QAC, is found to has a crucial bottleneck to solve a part of the optimization problem. Therefore we need to remove this problem while keeping its generality as a solver of the optimization problem. We again considered the application of JE to propose an alternative method, NQA. Although we do not assess its quantitative performance in the application to the optimization problem, our proposal gives a new paradigm to solve the optimization problem through the physical process like SA. We have to emphasize that QA was invented to solve the optimization problem in quantum computer. Therefore we must prepare the quantitative results to verify the precision and conditions in the actual experience on quantum computers. Along this line, we gave several results for the nonequilibrium behavior in the quantum system with gauge symmetry. These studies would be significant in the future development to realize the quantum computation.

Beyond the original version of SA, in order to find the desired solution as fast as possible, we must be away from the quasi-static procedure. The key point is to deal with nonequilibrium behavior. The further understanding of its peculiar behavior in statistical physics would be helpful to invent a genius and generic solver as PA and NQA.

Acknowledgement

The author thanks the fruitful discussions with Prof. Koji Hukushima of University of Tokyo and Prof. Hidetoshi Nishimori of Tokyo Institute of Technology. This work was partially supported by MEXT in Japan, Grant-in-Aid for Young Scientists (B) No.24740263.

© 2012 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Masayuki Ohzeki (August 29th 2012). Optimization by Use of Nature in Physics Beyond Classical Simulated Annealing, Simulated Annealing - Advances, Applications and Hybridizations, Marcos de Sales Guerra Tsuzuki, IntechOpen, DOI: 10.5772/50636. Available from:

chapter statistics

852total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Bayesian Recovery of Sinusoids with Simulated Annealing

By Dursun Üstündag and Mehmet Cevri

Related Book

First chapter

Mean Field Annealing Based Techniques for Resolving VLSI Automatic Design Problems

By G.R. Karimi and A. AziziVerki

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us