WO2018074006A1 - シミュレーション装置、コンピュータプログラム及びシミュレーション方法 - Google Patents

シミュレーション装置、コンピュータプログラム及びシミュレーション方法 Download PDF

Info

Publication number
WO2018074006A1
WO2018074006A1 PCT/JP2017/023378 JP2017023378W WO2018074006A1 WO 2018074006 A1 WO2018074006 A1 WO 2018074006A1 JP 2017023378 W JP2017023378 W JP 2017023378W WO 2018074006 A1 WO2018074006 A1 WO 2018074006A1
Authority
WO
WIPO (PCT)
Prior art keywords
magnetic field
magnetization
hamiltonian
function
probability distribution
Prior art date
Application number
PCT/JP2017/023378
Other languages
English (en)
French (fr)
Inventor
真之 大関
Original Assignee
国立大学法人京都大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立大学法人京都大学 filed Critical 国立大学法人京都大学
Priority to CA3041221A priority Critical patent/CA3041221A1/en
Priority to US16/341,024 priority patent/US20190235033A1/en
Publication of WO2018074006A1 publication Critical patent/WO2018074006A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/12—Measuring magnetic properties of articles or specimens of solids or fluids
    • G01R33/1284—Spin resolved measurements; Influencing spins during measurements, e.g. in spintronics devices
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N10/00—Quantum computing, i.e. information processing based on quantum-mechanical phenomena
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
    • G16Z99/00—Subject matter not provided for in other main groups of this subclass
    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03K—PULSE TECHNIQUE
    • H03K19/00—Logic circuits, i.e. having at least two inputs acting on one output; Inverting circuits
    • H03K19/02—Logic circuits, i.e. having at least two inputs acting on one output; Inverting circuits using specified components
    • H03K19/195—Logic circuits, i.e. having at least two inputs acting on one output; Inverting circuits using specified components using superconductive devices
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/044—Recurrent networks, e.g. Hopfield networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/047—Probabilistic or stochastic networks

Definitions

  • the present invention relates to a simulation apparatus, a computer program, and a simulation method.
  • optimization problem The problem of minimizing a discrete multivariable monovalent function (cost function) is called an optimization problem.
  • cost function The problem of minimizing a discrete multivariable monovalent function (cost function) is called an optimization problem.
  • Many important issues that require large-scale calculations such as pattern recognition, natural language processing, artificial intelligence, and machine learning can be formulated as optimization problems.
  • quantum annealing has attracted attention as an algorithm for solving optimization problems by skillfully utilizing the properties of quantum mechanics such as quantum fluctuations.
  • Quantum annealing expresses the cost function as an Ising model that is a function of binary variables, and formulates it as a problem to find the minimum value of the function. Quantum annealing is described in Non-Patent Document 1, for example.
  • Hamiltonian H0 expresses the cost function of the optimization problem, and H0 is selected so that the ground state of Hamiltonian H0 is the optimal solution.
  • ⁇ is a spin variable
  • the sigma of the x-direction component of the spin is an initial Hamiltonian representing a transverse magnetic field
  • the coefficient ⁇ is a parameter that controls the strength of quantum fluctuation.
  • the coefficient ⁇ is set to a very large value, and the coefficient ⁇ is decreased with the passage of time, and finally set to zero.
  • the state is explored by superimposing many states by large quantum fluctuations.
  • FIG. 15 is a schematic diagram showing an example of an energy gap of a quantum system.
  • the horizontal axis represents time, and the vertical axis represents energy.
  • ⁇ represents the energy gap between the ground state and the first excited state.
  • T represents the time required for quantum annealing, that is, the calculation time required to obtain the optimal solution.
  • N represents the number of spins.
  • the quantum phase transition first order phase transition
  • the calculation time T increases exponentially and becomes a very large value. For this reason, some quantum annealing has a very long calculation time. Even with ordinary computers, there are problems that require a very long time to solve the optimization problem, and applying quantum annealing to these problems similarly increases the calculation time. This is because of the existence of this quantum phase transfer.
  • Non-Patent Document 1 describes a technique for avoiding the first-order phase transition described above. That is, as shown in Equation (2), a square term (also referred to as antiferromagnetic XX interaction) of the sum of x-direction components of spin is added. In equation (2), gamma ⁇ is a coefficient.
  • FIG. 16 is a schematic diagram showing an example of a phase diagram.
  • the horizontal axis is the coefficient gamma ⁇
  • the vertical axis is the reciprocal of the coefficient ⁇ .
  • Symbol QP indicates a quantum paramagnetic phase
  • symbol F indicates a ferromagnetic phase.
  • the broken line indicates the boundary between QP and F.
  • the solid line in the horizontal direction indicates the first-order phase transition line.
  • the Hamiltonian represented by the formula (1) when the coefficient ⁇ is changed from a very large value to a small value as indicated by a broken line indicated by a symbol A, the problem of the first-order phase transition is encountered.
  • the Hamiltonian represented by the expression (2) the first-order phase transition can be avoided as indicated by the solid line indicated by the symbol B.
  • the antiferromagnetic XX interaction of Equation (2) has a term that is the square of the sum of the x-direction components of the spin, and the effect due to the x-direction component of the spin is an effect specific to quantum mechanics, Basically expressed in complex numbers.
  • the result may be a square of a complex number, that is, a negative value, resulting in a so-called negative sign problem in which the Boltzmann weight necessary for performing stochastic sampling is negative, It was impossible to simulate on a normal computer.
  • the present disclosure has been made in view of such circumstances, and an object thereof is to provide a simulation apparatus, a computer program, and a simulation method that can solve the negative sign problem while avoiding the primary phase transition.
  • the simulation apparatus represents a Hamiltonian of a system composed of a plurality of spins that can take two values as an initial Hamiltonian and a target Hamiltonian, and sets the initial Hamiltonian to a large value in an initial state, and changes the time.
  • a magnetization calculation unit that calculates the magnetization in the direction, an initial Hamiltonian calculation unit that calculates a magnetic field function including primary and secondary terms of magnetization calculated by the magnetization calculation unit as the initial Hamiltonian, and the magnetization calculation unit.
  • the difference between the calculated magnetization and the average of the sum of the components in the specified direction of the spin is a variable
  • a first probability distribution function calculation unit that calculates a probability distribution function for the initial Hamiltonian using an exponential operator including a multiplication term of the delta function and the magnetic field function, and the first probability distribution function calculation unit Based on the probability distribution obtained by calculation, a spin variable update unit that updates the spin variables of each of the plurality of spins, and whether the system is in an equilibrium state based on the spin variables updated by the spin variable update unit
  • a determination unit that determines whether or not, a first magnetization calculation unit that calculates the magnetization in the predetermined direction in the equilibrium state when the determination unit determines that the system is in an equilibrium state, and the first magnetization calculation
  • the simulation apparatus includes an exponential function type calculation including a derivative of the magnetic field function by integrating the delta function included in the probability distribution function calculated by the first probability distribution function calculation unit.
  • a second probability distribution function calculation unit that calculates a probability distribution function for the Hamiltonian of the system using a child, and the spin variable update unit is configured to calculate the probability distribution function calculated by the second probability distribution function calculation unit.
  • the spin variable is updated based on a probability distribution for the Hamiltonian.
  • the spin variable update unit determines that the magnetic field is not in a steady state by the magnetic field determination unit
  • the magnetic field function included in the probability distribution function for the Hamiltonian of the system The spin variable is updated based on a probability distribution with respect to a Hamiltonian of the system in which a derivative is updated based on the magnetization calculated by the first magnetization calculation unit.
  • a simulation apparatus is based on a setting unit that sets a plurality of magnetic field values in the predetermined direction in advance, the magnetic field value set in the setting unit, and an inverse function of a derivative of the magnetic field function.
  • a second magnetization calculation unit that calculates the magnetization in the predetermined direction, and the spin variable update unit includes a magnetic field set by the setting unit as a derivative of the magnetic field function included in a probability distribution function for a Hamiltonian of the system.
  • the spin variable is updated based on the probability distribution for the Hamiltonian of the system to which the value is assigned, and the determination unit determines whether the system is in an equilibrium state based on the spin variable updated by the spin variable update unit.
  • the first magnetization calculation unit calculates the magnetization in the predetermined direction in the equilibrium state
  • the physical quantity calculation unit When the magnetization calculated in the first magnetization calculator and the second magnetization calculator is equal, and calculates the physical quantity related to the system.
  • the computer program represents a Hamiltonian of a system composed of a plurality of spins that can take binary values as an initial Hamiltonian and a target Hamiltonian, and sets the initial Hamiltonian to a large value in an initial state
  • a process for calculating magnetization a process for calculating a magnetic field in the predetermined direction with respect to the plurality of spins based on the calculated magnetization and the magnetic field function, and a process for determining whether or not the calculated magnetic field is in a steady state.
  • a process for calculating a physical quantity related to the system is executed.
  • a Hamiltonian of a system composed of a plurality of spins that can take two values is represented by an initial Hamiltonian and a target Hamiltonian, and the initial Hamiltonian is set to a large value in an initial state, and the time variation is changed. Accordingly, a simulation method for simulating a physical quantity in an equilibrium state of the system so that the initial Hamiltonian is smaller than the target Hamiltonian, wherein an average of sums of predetermined direction components of the plurality of spins is calculated as the predetermined Hamiltonian.
  • a process for calculating the magnetization in the direction a process for calculating a magnetic field function including a first-order term and a second-order or higher-order term for the calculated magnetization as the initial Hamiltonian, an average of the sum of the calculated magnetization and a predetermined direction component of the spin, An exponential function including a multiplication term of the delta function with the difference between A process of calculating a probability distribution function for the initial Hamiltonian using a type operator, a process of updating the spin variables of each of the plurality of spins based on the probability distribution obtained by the calculation, and an updated spin variable
  • a process for determining whether the system is in an equilibrium state based on the above a process for calculating the magnetization in the predetermined direction in the equilibrium state when it is determined that the system is in an equilibrium state, Based on the magnetic field function, processing for calculating the magnetic field in the predetermined direction for the plurality of spins, processing for determining whether the calculated magnetic field is in a steady state, and determining that the magnetic field is in a steady state And
  • FIG. 1 is an explanatory diagram showing an example of the configuration of the simulation apparatus 100 according to the present embodiment.
  • the simulation apparatus 100 of the present embodiment can greatly expand the range of simulation by conventional quantum annealing, and can execute a simulation for solving an optimal solution using an Ising model.
  • the simulation apparatus 100 includes a control unit 10 that controls the entire apparatus, an input unit 11, a magnetization calculation unit 12, an initial Hamiltonian calculation unit 13, a first probability distribution function calculation unit 14, a second probability distribution function calculation unit 15, and a spin variable update.
  • Unit 16 equilibrium state determination unit 17, output unit 18, first magnetization calculation unit 19, magnetic field calculation unit 20, magnetic field determination unit 21, physical quantity calculation unit 22, second magnetization calculation unit 23, storage unit 24, and the like.
  • the input unit 11 receives input data (for example, the number of trotters, the number of spins, the value of the transverse magnetic field, etc.) for executing the simulation.
  • input data for example, the number of trotters, the number of spins, the value of the transverse magnetic field, etc.
  • the output unit 18 outputs output data (for example, energy E, magnetization m, etc.) as a result of the simulation.
  • the magnetization calculation unit 12 calculates the average of the sums of the plurality of spins in the predetermined direction as the magnetization in the predetermined direction.
  • the predetermined direction can be the x direction, which is the lateral direction, and the magnetization in the predetermined direction is the transverse direction.
  • Magnetization mx is the transverse direction.
  • the predetermined direction will be described as the horizontal direction.
  • the magnetization calculator 12 calculates the transverse magnetization mx as the sigma average of the x-direction component ⁇ x of the spin.
  • the simulation apparatus 100 is obtained by multiplying the function f having the variable of the sigma average of the x-direction component ⁇ x of the spin by the number N of spins as shown in the second term on the right side in the equation (3). Is formulated as an initial Hamiltonian.
  • the initial Hamiltonian acts as a quantum mechanical fluctuation that reverses the spin directed in the z direction.
  • the initial Hamiltonian computing unit 13 computes a magnetic field function including a transverse magnetization mx primary term and a second-order or higher term as an initial Hamiltonian.
  • the magnetic field function is represented by f (mx)
  • the coefficient ⁇ is a parameter for controlling the strength of quantum fluctuation
  • ⁇ is a predetermined coefficient.
  • the initial Hamiltonian can be represented by a magnetic field function.
  • the magnetic field function is not limited to the equation (4).
  • the above-mentioned “second or higher term” means that only a second order term, a higher order term of a third order term or more in addition to a second order term, a higher order term of a third order term or more without including a second order term. Means.
  • Equation (3) Hamiltonian H0, which is the first term on the right side, expresses the cost function of the optimization problem, and H0 is selected so that the ground state of Hamiltonian H0 is the optimal solution.
  • ⁇ is a spin variable and can take a value of ⁇ 1.
  • the coefficient ⁇ is set to a very large value, and the coefficient ⁇ is decreased with the passage of time, and finally set to zero.
  • the state is explored by superimposing many states by large quantum fluctuations.
  • the instantaneous ground state at each time is continuously traced and ⁇ gradually decreases, the relative weight of the Hamiltonian H0 becomes larger than the initial Hamiltonian, and finally the ground state of the Hamiltonian H0 is reached. In this state, the solution to the optimization problem is obtained, and the required physical quantity can be calculated.
  • the first probability distribution function calculation unit 14 uses an exponential function type operator including a multiplication term of a delta function having a variable between the transverse magnetization mx and the average of the sum of the transverse components ⁇ x of the spin and a magnetic field function. To calculate the probability distribution function for the initial Hamiltonian. Specifically, the probability distribution function calculated by the first probability distribution function calculation unit 14 can be expressed as in Expression (5).
  • Equation (5) shows the probability distribution when Suzuki Trotter decomposition is performed.
  • ⁇ is a coefficient proportional to the reciprocal of absolute temperature, and ⁇ indicates the number of trotter.
  • Equation (5) by rewriting the effect on the magnetic field function f (mx), any quantum fluctuations including f (mx) (the problem of the initial Hamiltonian represented by the second term on the right side of Equation (3)) ) Can be changed to a simple one with a transverse magnetic field.
  • the transverse magnetization mx When equal to the average of the sum, the delta function is 1, and the multiplication term of the delta function and the magnetic field function f (mx) can be replaced by the magnetic field function f (mx) as a result, and the exponential function of the probability distribution function Since the second-order or higher-order term (including the XX interaction) of the sum of the x-direction components ( ⁇ x) of the spin can be removed from the type operator, the negative sign problem is eliminated.
  • Equation (6) represents the Suzuki Trotter decomposition formula.
  • a and B are operators, and L is the number of trotters.
  • the number of trotters is also represented by ⁇ .
  • a Hamiltonian of a quantum system can be defined by a sum of local Hamiltonians that represent interactions between local components. Since local Hamiltonians are not commutative with each other, the size of the matrix representation of the Hamiltonian of the quantum system increases, and the calculation cost increases. Therefore, by using the Suzuki Trotter decomposition formula shown in Equation (6), the exponent operator can be decomposed into a product of local Hamiltonian exponent operators having a small matrix representation size.
  • the second probability distribution function calculation unit 15 displays the delta function included in the probability distribution function calculated by the first probability distribution function calculation unit 14 as an integral display, and uses an exponential function type operator including a derivative of the magnetic field function. Compute the probability distribution function for the Hamiltonian of the system. Specifically, the probability distribution function calculated by the second probability distribution function calculation unit 15 can be expressed as Expression (7).
  • Equation (7) ⁇ is the number of trotters, and m tilde x can be represented by a derivative f ′ (mx) of mx of the magnetic field function f (mx) as in equation (8).
  • m tilde x is introduced by performing integral display (Fourier integral display) of the delta function of Expression (5).
  • Equation (7) the Hamiltonian problem having the x-direction component ⁇ x of the spin can be treated as a problem of the transverse magnetic field, and by executing the Suzuki Trotter decomposition, the term related to the x-direction component ⁇ x of the spin is expressed by Equation (9).
  • the coefficient B can be expressed by Expression (10).
  • Expression (9) is also referred to as a probability distribution function calculated by the second probability distribution function calculation unit 15.
  • Z (mx) is a coefficient for normalization.
  • FIG. 2 is a schematic diagram showing an example of Suzuki Trotter decomposition.
  • the horizontal axis represents a spin variable on a one-dimensional lattice and indicates a so-called real space direction.
  • the vertical axis is the direction (trotter direction) introduced by Suzuki Trotter decomposition, and state variables are arranged on two-dimensional lattice points.
  • the quantum model is transformed into a classical model having a state space with one more dimension by Suzuki Trotter decomposition.
  • Equation (9) can be derived by using the saddle point method that focuses only on the place where the maximum value of mx is taken.
  • the expression representing the saddle point condition (becomes a maximum point) is Expression (8), and m tilde x disappears in Expression (9).
  • Equation (9) shows how the spin variable ⁇ behaves when the transverse magnetization mx is determined.
  • the spin variable update unit 16 updates the spin variable of each of the plurality of spins based on the probability distribution obtained by the calculation by the first probability distribution function calculation unit 14. More specifically, the spin variable update unit 16 updates the spin variable based on the probability distribution for the Hamiltonian of the system obtained by the calculation by the second probability distribution function calculation unit 15.
  • the spin variable can be updated, for example, by selecting a spin variable and satisfying a detailed balance such as a heat bath method or a metropolis method.
  • FIG. 3 is a schematic diagram showing an example of how the spin variable is updated.
  • the left figure shows the current state of the spin variable, and the right figure assumes that the spin variable of the lattice i is selected as an arbitrary spin variable. For convenience, it is assumed that there are four spin variables around the spin variable of the lattice i, and only the spin variable of the lattice i can be changed (1 spin flip).
  • the current probability distribution Pp is calculated by substituting the spin variable in the state of the left figure into the equation (9).
  • the next new probability distribution Pn is calculated by substituting the spin variable in the state where the spin variable of the selected lattice i is changed into the equation (9).
  • the equilibrium state determination unit 17 has a function as a determination unit, and determines whether or not the system is in an equilibrium state based on the updated spin variable. To determine whether or not the system has reached an equilibrium state, for example, the system energy E and magnetization m (m-ment such as the square of the magnetization, the fourth power, etc.) are calculated (measured). It can be determined that an equilibrium state has been reached.
  • the energy E of the system is calculated by equation (1).
  • the magnetization m can be obtained by calculating the sum for all spins and dividing the sum by the number of spins.
  • the first magnetization calculation unit 19 calculates the transverse magnetization mx in the equilibrium state.
  • the transverse magnetization mx in the equilibrium state can be obtained by the time average value of the amount calculated by the equation (11).
  • i represents the spin location (lattice point)
  • t represents the number of trotters
  • ⁇ represents the total number of trotters
  • N represents the total number of spins.
  • the magnetic field calculation unit 20 calculates a transverse magnetic field for a plurality of spins based on the transverse magnetization mx and the magnetic field function f (mx) calculated by the first magnetization calculation unit 19.
  • the transverse magnetic field (m tilde x) can be calculated by substituting the transverse magnetization mx into the equation (8) obtained by differentiating the magnetic field function f (mx) by the transverse magnetization mx.
  • the magnetic field determination unit 21 determines whether or not the transverse magnetic field (m tilde x) calculated by the magnetic field calculation unit 20 is in a steady state. Whether the transverse magnetic field is in a steady state or not is determined by changing the value of the transverse magnetic field adaptively according to the value of the transverse magnetization mx in the equilibrium state, and determining that the transverse magnetic field is in a steady state when the transverse magnetic field does not change. Can do.
  • the physical quantity calculation unit 22 calculates a physical quantity related to the system when the magnetic field determination unit 21 determines that the transverse magnetic field is in a steady state.
  • the physical quantities related to the system are, for example, energy E and magnetization m. While continuing the simulation, physical quantities such as energy E and magnetization m are repeatedly calculated based on the spin variables, and when a certain amount of time has elapsed, the time average of the physical quantities is calculated and used as the result of the physical quantity. The time average can be calculated with an arbitrary accuracy, and if the time is lengthened, the accuracy can be increased accordingly.
  • the storage unit 24 stores data necessary for simulation, input data, processing results obtained during the simulation, output data, and the like.
  • the spin variable update unit 16 calculates the derivative f ′ (mx) of the magnetic field function f (mx) included in the probability distribution function for the system Hamiltonian.
  • the spin variable is updated based on the probability distribution for the Hamiltonian of the updated system, updated based on the transverse magnetization calculated by the first magnetization calculator 19.
  • the coefficient B in the equation (9) includes the derivative f ′ (mx) of the magnetic field function f (mx) as shown by the equation (10). That is, when f ′ (mx) is updated by mx, the coefficient B changes, and as a result, the probability distribution of the system calculated by Expression (9) changes.
  • FIG. 4 is an explanatory diagram showing an example of the relationship between the transverse magnetization mx and the transverse magnetic field m tilde x.
  • m tilde x is a derivative f ′ (mx) of the magnetic field function f (mx).
  • the derivative function f ′ (mx) is a linear function of the transverse magnetization mx. That is, the transverse magnetic field (m tilde x) can be changed by changing the transverse magnetization mx. In the conventional quantum annealing, the transverse magnetic field was constant.
  • the solution can be reached by calculating the probability distribution for the Hamiltonian of the system according to the transverse magnetization and performing the process of updating the spin variable again.
  • FIG. 5 is a flowchart showing an example of the processing procedure of the adaptive quantum Monte Carlo method performed by the simulation apparatus 100 according to the present embodiment.
  • “Adaptive” in the adaptive quantum Monte Carlo method is a term used to distinguish it from the quantum Monte Carlo method, which is a probabilistic method for realizing the conventional quantum annealing method.
  • the negative sign problem can be avoided while avoiding, and this represents a simulation method that can be executed in a normal computer.
  • the processing subject will be described as the control unit 10.
  • Control unit 10 sets the number of trotter and the number of spins (S11).
  • the number of trotters depends on the performance of the simulator (computer), but can be 128, for example.
  • the number N of spins can be set to an arbitrary size. If the number of trotter and the number of spins are increased, the calculation accuracy can be improved.
  • the control unit 10 sets an initial value for the sigma of the spin variable (S12). If an initial value is set for the sigma of the spin variable, the transverse magnetization mx can be calculated, so the value of the transverse magnetic field (m tilde x, ie, f '(mx)) can be set to the initial value.
  • the control unit 10 initializes the spin variable with a random number (S13), selects the spin variable, and updates it with the heat bath method or the metropolis method (that satisfies the detailed balance) (S14). This causes the Hamiltonian of the system to converge toward the ground state of the target Hamiltonian H0.
  • Control unit 10 determines whether or not the system is in an equilibrium state (S15). The determination as to whether or not an equilibrium state has been reached can be made based on whether or not the value of energy E, the magnetization m, etc. fluctuate.
  • control unit 10 continues the processing from step S14.
  • the control unit 10 changes the value of the transverse magnetic field (m tilde x, that is, f ′ (mx)) according to the value of the transverse magnetization mx in the equilibrium state ( S16).
  • the control unit 10 determines whether the transverse magnetic field changes, that is, whether the transverse magnetic field is in a steady state (S17). When the transverse magnetic field changes (YES in S17), the control unit 10 continues the processing from step S14 onward, assuming that another transverse magnetic field or transverse magnetization that is not a solution is obtained.
  • control unit 10 starts measuring (calculating) the physical quantity assuming that the solution is obtained (S18), and when the required time has elapsed, the time average of the measured quantity Is obtained as the result of the physical quantity, and the process is terminated.
  • FIG. 6 is an explanatory diagram showing a first example of a simulation result by the adaptive quantum Monte Carlo method of the present embodiment.
  • the horizontal axis indicates the horizontal magnetic field ⁇ (Gamma), and the vertical axis indicates the energy E.
  • N indicates the number of spins during simulation
  • the simulation result almost coincides with the true solution, and in particular, it can be seen that the closer to the true solution by increasing the number of spins N.
  • FIG. 7 is an explanatory diagram showing a second example of the simulation result by the adaptive quantum Monte Carlo method of the present embodiment
  • FIG. 8 is an explanation showing a third example of the simulation result by the adaptive quantum Monte Carlo method of the present embodiment.
  • FIG. 7 the vertical axis indicates the magnetization m.
  • the vertical axis represents the transverse magnetization mx. 7 and 8, it can be seen that the simulation result is closer to the true solution by increasing the spin number N.
  • the quantum Monte Carlo method based on theta analysis is an ordinary quantum Monte Carlo method.
  • the control unit 10 has a function as a setting unit, and sets a plurality of transverse magnetic field values in advance. That is, a plurality of transverse magnetic field values are prepared.
  • the spin variable update unit 16 updates the spin variable based on the probability distribution for the system Hamiltonian in which the value of the transverse magnetic field set by the control unit 10 is assigned to the derivative of the magnetic field function included in the probability distribution function for the system Hamiltonian. To do.
  • the equilibrium state determination unit 17 determines whether or not the system is in an equilibrium state based on the spin variable updated by the spin variable update unit 16.
  • the first magnetization calculator 19 calculates the transverse magnetization in the equilibrium state when the system is in an equilibrium state. That is, the transverse magnetization is calculated by executing the quantum Monte Carlo method using a previously facilitated transverse magnetic field. Thereby, the relationship between a transverse magnetic field and transverse magnetization can be plotted.
  • the physical quantity calculator 22 calculates a physical quantity related to the system when the transverse magnetization calculated by the first magnetization calculator 19 and the second magnetization calculator 23 is equal. That is, when the transverse magnetization calculated based on the inverse function of the derivative of the magnetic field function is equal to the transverse magnetization calculated by executing the quantum Monte Carlo method, by obtaining the transverse magnetization and the corresponding transverse magnetic field, A solution can also be obtained by a data analysis approach.
  • FIG. 9 is a flowchart showing an example of a processing procedure of the quantum Monte Carlo method by data analysis performed by the simulation apparatus 100 of the present embodiment.
  • the control unit 10 sets a plurality of values of the transverse magnetic field (m tilde x, i.e., f '(mx)) (S31), and executes the quantum Monte Carlo method using each of the set values of the transverse magnetic field (S32). .
  • m tilde x i.e., f '(mx)
  • the control unit 10 associates the value of the transverse magnetization obtained by executing the quantum Monte Carlo method with the transverse magnetic field when the transverse magnetization is obtained, and the relationship between the value of the transverse magnetic field and the value of the transverse magnetization. Is plotted (S33). Here, plotting does not have to be actually drawn on a chart, but may be in a form in which the correspondence can be understood.
  • the control unit 10 calculates transverse magnetization based on the inverse function of the derivative f ′ (mx) of the magnetic field function f (mx) (S34). Specifically, the transverse magnetization can be calculated by substituting each set value of the transverse magnetic field into an inverse function.
  • the control unit 10 specifies the value of the transverse magnetization and the value of the transverse magnetic field when the transverse magnetization calculated by the inverse function matches the transverse magnetization on the plot (S35). As a result, a true solution can be obtained, and a physical quantity result can be obtained in the same manner as in the adaptive quantum Monte Carlo method.
  • the control unit 10 ends the process.
  • the value of the transverse magnetization mx and the value of the transverse magnetic field m tilde x at the point where the curve P1 and the straight line P2 intersect are the solutions.
  • FIG. 11 is an explanatory diagram showing a first example of a simulation result by a quantum Monte Carlo method based on data analysis of the present embodiment.
  • the horizontal axis represents the transverse magnetic field ⁇ (Gamma), and the vertical axis represents the energy E.
  • N indicates the number of spins during simulation
  • the chart corresponding to each N is a set of points where the curve P1 and the straight line P2 intersect.
  • the simulation result almost coincides with the true solution, and it can be seen that the simulation solution approaches the true solution by increasing the number of spins N in particular.
  • FIG. 12 is an explanatory view showing a second example of the simulation result by the quantum Monte Carlo method based on the data analysis of the present embodiment
  • FIG. 13 is a third example of the simulation result by the quantum Monte Carlo method by the data analysis of the present embodiment. It is explanatory drawing shown.
  • the vertical axis indicates the magnetization m.
  • the vertical axis represents the transverse magnetization mx. 12 and 13, it can be seen that the simulation result is closer to a true solution by increasing the spin number N.
  • the simulation apparatus 100 can be realized using a computer including a CPU (processor), a RAM (memory), and the like. That is, as shown in FIG. 5 and FIG. 9, a computer program that defines the procedure of each process is loaded into a RAM (memory) provided in the computer, and the computer program is executed by a CPU (processor). 100 can be realized.
  • a computer program that defines the procedure of each process is loaded into a RAM (memory) provided in the computer, and the computer program is executed by a CPU (processor). 100 can be realized.
  • FIG. 14 is an explanatory diagram showing another example of the configuration of the simulation apparatus of the present embodiment.
  • reference numeral 300 denotes a normal computer.
  • the computer 300 includes a control unit 30, an input unit 40, an output unit 50, an external I / F (interface) unit 60, and the like.
  • the control unit 30 includes a CPU 31, a ROM 32, a RAM 33, an I / F (interface) 34, and the like.
  • the input unit 40 acquires input data for simulation.
  • the output unit 50 outputs output data that is a simulation result.
  • the I / F 34 has an interface function between the control unit 30 and each of the input unit 40, the output unit 50, and the external I / F unit 60.
  • the external I / F unit 60 can read a computer program from a recording medium M (for example, a medium such as a DVD) on which the computer program is recorded.
  • a recording medium M for example, a medium such as a DVD
  • the computer program recorded on the recording medium M is not limited to the one recorded on a portable medium, and may be a computer program transmitted through the Internet or another communication line. Can be included.
  • the computer includes a computer system including a single computer having a plurality of processors or a plurality of computers connected via a communication network.
  • the present invention is not limited to a limited model, and is not limited to a wide range of models. Simulation can be executed on a computer, the application range of the quantum Monte Carlo method can be expanded, and the scope of material search or design simulation can be expanded. In addition, the simulation method of the present embodiment can be used in technical fields that require large-scale calculations such as the production site of quantum computers, artificial intelligence, and machine learning.
  • the simulation apparatus represents a Hamiltonian of a system composed of a plurality of spins that can take binary values as an initial Hamiltonian and a target Hamiltonian, and sets the initial Hamiltonian to a large value in an initial state, and with time change
  • a magnetization calculation unit that calculates the magnetization, an initial Hamiltonian calculation unit that calculates, as the initial Hamiltonian, a magnetic field function including primary and secondary terms of magnetization calculated by the magnetization calculation unit, and a calculation performed by the magnetization calculation unit
  • the variable is the difference between the magnetization and the average of the sums of the spin components in the specified direction
  • a first probability distribution function computing unit that computes a probability distribution function for the initial Hamiltonian using an exponential function operator including a multiplication term of a ruther function and the magnetic field function; and computing by the first probability distribution function computing unit
  • a spin variable update unit that updates the spin variables of each of the plurality of spins based on the probability distribution obtained as a result, and whether the system is in an equilibrium state based on the spin variables updated by the spin variable update unit
  • a determination unit that determines whether the system is in an equilibrium state, and a first magnetization calculation unit that calculates the magnetization in the predetermined direction in the equilibrium state, and the first magnetization calculation unit
  • the computer program according to the present embodiment represents a Hamiltonian of a system composed of a plurality of spins that can take binary values as an initial Hamiltonian and a target Hamiltonian, and sets the initial Hamiltonian to a large value in an initial state, and changes with time.
  • a process of calculating the average as the magnetization in the predetermined direction a process of calculating a magnetic field function including a first-order term and a second-order or higher term of the calculated magnetization as the initial Hamiltonian, and the calculated magnetization and the predetermined direction component of the spin A delta function whose variable is the difference from the average of the sum, and A process for calculating a probability distribution function for the initial Hamiltonian using an exponential operator including a multiplication term with a field function, and a spin variable of each of the plurality of spins based on the probability distribution obtained by the calculation
  • a process for updating a process for determining whether or not the system is in an equilibrium state based on the updated spin variable, and a magnetization in the predetermined direction in the equilibrium state when it is determined that the system is in an equilibrium state
  • a process for calculating the magnetic field in the predetermined direction for the plurality of spins based on the calculated magnetization and the magnetic field function a process for determining whether the calculated magnetic field is in a steady state
  • the computer represents a Hamiltonian of a system composed of a plurality of spins that can take binary values as an initial Hamiltonian and a target Hamiltonian, and the initial Hamiltonian has a large value in an initial state.
  • a Hamiltonian of a system composed of a plurality of spins that can take binary values is represented by an initial Hamiltonian and a target Hamiltonian, and the initial Hamiltonian is set to a large value in an initial state, and with time change
  • the magnetization calculation unit calculates the average of the sums of the plurality of spins in the predetermined direction as the magnetization in the predetermined direction.
  • the predetermined direction can be the x direction, which is the lateral direction, and the magnetization in the predetermined direction is the transverse direction.
  • Magnetization mx That is, the transverse magnetization mx is calculated as the sigma average of the x-direction component of the spin.
  • the initial Hamiltonian calculation unit calculates a magnetic field function including the calculated first and second-order terms of magnetization as an initial Hamiltonian.
  • the magnetic field function is represented by f (mx)
  • the initial Hamiltonian can be represented by a magnetic field function.
  • the first probability distribution function calculation unit is initialized using an exponential operator including a multiplication term of a delta function having a variable between a calculated magnetization and an average of sums of components in a predetermined direction of spin and a magnetic field function. Compute the probability distribution function for the Hamiltonian.
  • the delta function is 1 Since the multiplication term of the delta function and the magnetic field function is replaced with the magnetic field function as a result, the higher-order term of the second or higher order of the sum of the predetermined direction components ( ⁇ x) of the spin from the exponential function operator of the probability distribution function is Since it can be eliminated, the minus sign problem is eliminated.
  • the spin variable update unit updates the spin variable of each of the plurality of spins based on the probability distribution obtained by the calculation. For example, a spin variable can be selected and updated with a detailed balance such as a hot bath method or a metropolis method.
  • the determination unit determines whether or not the system is in an equilibrium state based on the updated spin variable. To determine whether the system is in an equilibrium state, for example, calculate (measure) the energy and magnetization of the system (m-ment such as the square of the magnetization, the fourth power, etc.), and if there is no fluctuation, the system is in an equilibrium state. It can be determined that
  • the first magnetization calculator calculates the magnetization (transverse magnetization mx) in a predetermined direction in the equilibrium state when the system is in an equilibrium state.
  • the magnetic field calculation unit calculates a magnetic field (transverse magnetic field) in a predetermined direction for a plurality of spins based on the calculated magnetization (transverse magnetization) and the magnetic field function.
  • the transverse magnetic field can be calculated by substituting transverse magnetization into an expression obtained by differentiating the magnetic field function by transverse magnetization.
  • the magnetic field determination unit determines whether or not the calculated magnetic field (transverse magnetic field) is in a steady state. Whether the transverse magnetic field is in a steady state can be determined as a steady state when the transverse magnetic field does not change by adaptively changing the value of the transverse magnetic field according to the value of transverse magnetization in the equilibrium state. it can.
  • the physical quantity calculation unit calculates a physical quantity related to the system when it is determined that the transverse magnetic field is in a steady state.
  • the physical quantities related to the system are, for example, energy and magnetization. While continuing the simulation, the physical quantities such as energy and magnetization are repeatedly calculated based on the spin variables, and the time average of the physical quantities is calculated when a certain amount of time has passed, and the result is the physical quantity.
  • the delta function included in the probability distribution function calculated by the first probability distribution function calculation unit is displayed as an integral display, and an exponential function type operator including the derivative of the magnetic field function is displayed.
  • a second probability distribution function calculation unit that calculates a probability distribution function for the Hamiltonian of the system using the second probability distribution function calculation unit, and the spin variable update unit for the Hamiltonian of the system obtained by the calculation by the second probability distribution function calculation unit The spin variable is updated based on the probability distribution.
  • the second probability distribution function calculation unit displays the delta function included in the probability distribution function calculated by the first probability distribution function calculation unit as an integral display, and uses an exponential type operator including a derivative of the magnetic field function. Compute the probability distribution function for the Hamiltonian.
  • the delta function integral the derivative of the magnetic field function can be introduced, and by performing the Suzuki Trotter decomposition, the term relating to the x-direction component ( ⁇ x) of the spin can be replaced with the interaction between the trotter. It is possible to perform numerical calculation.
  • the spin variable update unit updates the spin variable based on the probability distribution for the Hamiltonian of the system obtained by the calculation by the second probability distribution function calculation unit. Thereby, simulation can be performed using a normal computer.
  • the spin variable update unit determines that the magnetic field is not in a steady state by the magnetic field determination unit
  • the derivative of the magnetic field function included in the probability distribution function for the Hamiltonian of the system Is updated based on the probability distribution for the Hamiltonian of the system updated based on the magnetization calculated by the first magnetization calculator.
  • the spin variable update unit determines that the transverse magnetic field is not in a steady state
  • the derivative of the magnetic field function included in the probability distribution function for the Hamiltonian of the system is updated based on the transverse magnetization calculated by the first magnetization calculation unit.
  • the spin variable is updated based on the probability distribution for the Hamiltonian.
  • the solution can be reached by calculating the probability distribution for the Hamiltonian of the system according to the transverse magnetization and performing the process of updating the spin variable again.
  • the simulation apparatus includes a setting unit that sets a plurality of magnetic field values in the predetermined direction in advance, the magnetic field value set by the setting unit, and an inverse function of a derivative of the magnetic field function.
  • a second magnetization calculation unit that calculates a magnetization in the direction, and the spin variable update unit sets a magnetic field value set by the setting unit to a derivative of the magnetic field function included in a probability distribution function for a Hamiltonian of the system
  • the spin variable is updated based on the assigned probability distribution for the Hamiltonian of the system, and the determination unit determines whether the system is in an equilibrium state based on the spin variable updated by the spin variable update unit.
  • the first magnetization calculation unit calculates the magnetization in the predetermined direction in the equilibrium state
  • the physical quantity calculation unit When the magnetization calculated in the first magnetization calculator and the second magnetization calculator is equal, and calculates the physical quantity related to the system.
  • the setting unit sets a plurality of values of a magnetic field (transverse magnetic field) in a predetermined direction in advance. That is, a plurality of transverse magnetic field values are prepared.
  • the second magnetization calculator calculates the magnetization (transverse magnetization) in a predetermined direction based on the value of the transverse magnetic field and the inverse function of the derivative of the magnetic field function. By calculating the transverse magnetization for a plurality of transverse magnetic fields, the relationship between the transverse magnetic field and the transverse magnetization can be plotted.
  • the spin variable update unit updates the spin variable based on the probability distribution for the system Hamiltonian in which the magnetic field value set by the setting unit is assigned to the derivative of the magnetic field function included in the probability distribution function for the system Hamiltonian.
  • the determination unit determines whether or not the system is in an equilibrium state based on the spin variable updated by the spin variable update unit.
  • the first magnetization calculation unit calculates the transverse magnetization in the equilibrium state. That is, the transverse magnetization is calculated by executing the quantum Monte Carlo method using a previously facilitated transverse magnetic field. Thereby, the relationship between a transverse magnetic field and transverse magnetization can be plotted.
  • the physical quantity calculation unit calculates a physical quantity related to the system when the magnetizations calculated by the first magnetization calculation unit and the second magnetization calculation unit are equal. That is, when the transverse magnetization calculated based on the inverse function of the derivative of the magnetic field function is equal to the transverse magnetization calculated by executing the quantum Monte Carlo method, by obtaining the transverse magnetization and the corresponding transverse magnetic field, A solution can also be obtained by a data analysis approach.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Computer Hardware Design (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computational Linguistics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Mram Or Spin Memory Techniques (AREA)

Abstract

一次相転移を回避しつつ負符号問題を解決することができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供する。　複数のスピンの所定方向成分の和の平均を所定方向の磁化として演算する磁化演算部と、磁化の一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数と磁場関数との乗算項を含む、初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、スピン変数を更新するスピン変数更新部と、系が平衡状態になったか否かを判定する判定部と、系の平衡状態での所定方向の磁化を算出する第１磁化算出部と、所定方向の磁場を算出する磁場算出部と、磁場が定常状態であるか否かを判定する磁場判定部と、系に関する物理量を算出する物理量算出部とを備える。

Description

シミュレーション装置、コンピュータプログラム及びシミュレーション方法
　本発明は、シミュレーション装置、コンピュータプログラム及びシミュレーション方法に関する。
　離散的な多変数の１価関数（コスト関数）を最小化する問題を最適化問題という。パターン認識、自然言語処理、人工知能、機械学習をはじめとする大規模な計算を必要とする多くの重要な課題が最適化問題として定式化できる。また、近年、量子揺らぎなどの量子力学の性質を巧みに利用して最適化問題を解くアルゴリズムとしての量子アニーリングが注目されている。
　量子アニーリングでは、コスト関数を二値変数の関数であるイジング（Ising）モデルとして表現し、その関数の最小値を求める問題として定式化する。量子アニーリングは、例えば、非特許文献１に記載されている。
　量子アニーリングでは、式（１）で示すように、ｚ方向に向いたスピン（ベクトル表記法で示す太文字σで表すように多くのスピン）を反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する。この量子揺らぎの効果によって、より最適な解の探索を行う。
Figure JPOXMLDOC01-appb-I000001
　ここで、ハミルトニアンＨ０は、最適化問題のコスト関数を表現したものであり、ハミルトニアンＨ０の基底状態が最適解となるようにＨ０を選択する。σはスピン変数であり、スピンのｘ方向成分のシグマは、横磁場を表す初期ハミルトニアンであり、係数Γは、量子揺らぎの強さを制御するパラメータである。初期状態（時刻ｔ＝０）では、係数Γを非常に大きな値とし、時間の経過とともに係数Γを小さくし、最終的には０にする。最初は大きな量子揺らぎによって数多くの状態の重ね合わせを実現して状態探査をする。各時刻における瞬間的な基底状態を連続的にたどり、次第にΓが小さくなると、初期ハミルトニアンに比べてハミルトニアンＨ０の相対的な重みが大きくなり、最終的には、ハミルトニアンＨ０の基底状態に到達する。これは、最適化問題の解が得られ、所要の物理量が算出することができたということである。
　図１５は量子系のエネルギーギャップの一例を示す模式図である。図１５において、横軸は時間を示し、縦軸はエネルギーを示す。量子系のエネルギー準位を時系列で追いかけた場合、基底状態と励起状態のエネルギー準位が接近する場合がある。Δは基底状態と第一励起状態とのエネルギーギャップを表す。Ｔは量子アニーリングに要する時間、すなわち最適解を求めるまでに要する計算時間を示す。また、Ｎはスピンの数を示す。図１５に示すようなエネルギーギャップΔでは、量子相転移（一次相転移）が起こり、計算時間Ｔは、指数関数的に増大して非常に大きな値となる。このため、量子アニーリングでも非常に計算時間が長いものも存在する。通常のコンピュータでも最適化問題を解くために非常に長い時間を要する問題があり、それらの問題に量子アニーリングを適用すると、同様に計算時間が長くなる。その背景には、この量子相移転の存在があるためである。
　非特許文献１には、上述の一次相転移を回避する手法が記載されている。すなわち、式（２）に示すように、スピンのｘ方向成分の和の２乗の項（反強磁性ＸＸ相互作用とも称する）を追加する。式（２）において、ガンマγは係数である。
Figure JPOXMLDOC01-appb-I000002
　図１６は相図の一例を示す模式図である。図１６において、横軸は係数ガンマγであり、縦軸は係数Γの逆数である。符号ＱＰは量子常磁性相を示し、符号Ｆは強磁性相を示す。破線はＱＰとＦとの境界を示す。また、図１６中、横方向の実線は、一次相転移の線を示す。式（１）で示すハミルトニアンの場合には、符号Ａで示す破線のように、係数Γを非常に大きな値から小さい値にする場合、一次相転移の問題に遭遇する。一方、式（２）で示すハミルトニアンの場合には、符号Ｂで示す実線のように、一次相転移を回避することができる。
Yuya Seki, Hidetoshi Nishimori, "Quantum annealing with antiferromagnetic fluctuations", Phys. Rev. E 85, 051112(2012)
　しかし、式（２）の反強磁性ＸＸ相互作用は、スピンのｘ方向成分の和の２乗の項を有しており、スピンのｘ方向成分による効果は、量子力学特有の効果であり、基本的には複素数で表現される。ＸＸ相互作用を利用した場合、その結果に複素数の２乗、すなわち負の値を生じる場合があり、確率的なサンプリングを行う際に必要なボルツマン重みが負となる、いわゆる負符号問題が生じ、通常のコンピュータ上でのシミュレーションが不可能とされていた。
　本開示は斯かる事情に鑑みてなされたものであり、一次相転移を回避しつつ負符号問題を解決することができるシミュレーション装置、コンピュータプログラム及びシミュレーション方法を提供することを目的とする。
　本開示の実施の形態に係るシミュレーション装置は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部とを備えることを特徴とする。
　本開示の実施の形態に係るシミュレーション装置は、前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、前記スピン変数更新部は、前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
　本開示の実施の形態に係るシミュレーション装置は、前記スピン変数更新部は、前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
　本開示の実施の形態に係るシミュレーション装置は、予め前記所定方向の磁場の値を複数設定する設定部と、該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部とを備え、前記スピン変数更新部は、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、前記判定部は、前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、前記第１磁化算出部は、前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、前記物理量算出部は、前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする。
　本開示の実施の形態に係るコンピュータプログラムは、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させることを特徴とする。
　本開示の実施の形態に係るシミュレーション方法は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを含むことを特徴とする。
　本開示によれば、一次相転移を回避しつつ負符号問題を解決することができる。
本実施の形態のシミュレーション装置の構成の一例を示す説明図である。 鈴木トロッター分解の一例を示す模式図である。 スピン変数の更新の様子の一例を示す模式図である。 横磁化ｍｘと横磁場ｍチルダｘとの関係の一例を示す説明図である。 本実施の形態のシミュレーション装置による適応的量子モンテカルロ法の処理手順の一例を示すフローチャートである。 本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。 本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図である。 本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。 本実施の形態のシミュレーション装置が行うデータ解析による量子モンテカルロ法の処理手順の一例を示すフローチャートである。 本実施のデータ解析による量子モンテカルロ法の概念を示す説明図である。 本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。 本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図である。 本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。 本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。 量子系のエネルギーギャップの一例を示す模式図である。 相図の一例を示す模式図である。
　以下、本実施の形態を図面に基づいて説明する。図１は本実施の形態のシミュレーション装置１００の構成の一例を示す説明図である。本実施の形態のシミュレーション装置１００は、従来の量子アニーリングによるシミュレーションの範囲を大幅に広げることができ、イジング（Ising）モデルを用いて最適解を解くためのシミュレーションを実行することができる。
　シミュレーション装置１００は、装置全体を制御する制御部１０、入力部１１、磁化演算部１２、初期ハミルトニアン演算部１３、第１確率分布関数演算部１４、第２確率分布関数演算部１５、スピン変数更新部１６、平衡状態判定部１７、出力部１８、第１磁化算出部１９、磁場算出部２０、磁場判定部２１、物理量算出部２２、第２磁化算出部２３、記憶部２４などを備える。
　入力部１１は、シミュレーションを実行するための入力データ（例えば、トロッター数、スピンの数、横磁場の値など）を受け付ける。
　出力部１８は、シミュレーションの結果である出力データ（例えば、エネルギーＥ、磁化ｍなど）を出力する。
　磁化演算部１２は、複数のスピンの所定方向成分の和の平均を所定方向の磁化として演算する。ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する場合、所定方向は、横方向であるｘ方向とすることができ、所定方向の磁化は、横方向の磁化ｍｘとすることができる。以下では、所定方向を横方向として説明する。磁化演算部１２は、横磁化ｍｘを、スピンのｘ方向成分σｘのシグマの平均として演算する。
　本実施の形態のシミュレーション装置１００は、式（３）において、右辺第２項で示すように、スピンのｘ方向成分σｘのシグマの平均を変数とする関数ｆにスピンの数Ｎを乗算したものを、初期ハミルトニアンとして定式化する。初期ハミルトニンアンは、ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして作用する。
Figure JPOXMLDOC01-appb-I000003
　具体的には、初期ハミルトニアン演算部１３は、横磁化ｍｘ一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する。式（４）に示すように、磁場関数をｆ（ｍｘ）で表すと、磁場関数ｆ（ｍｘ）は、例えば、ｆ（ｍｘ）＝Γ・ｍｘ＋（γ／２）・（ｍｘ）2 と表すことができる。ここで、係数Γは、量子揺らぎの強さを制御するパラメータであり、γは所定の係数である。初期ハミルトニアンは、磁場関数で表すことができる。磁場関数に横磁化ｍｘの２乗の項を含めることにより、一次相転移の問題を回避することが可能となる。なお、磁場関数は、式（４）に限定されるものではない。なお、上述の「二次以上の項」とは、二次項のみの場合、二次項に加えて三次項以上の高次項を含む場合、二次項は含まずに三次項以上の高次項を含む場合を意味する。
　式（３）において、右辺の第１項である、ハミルトニアンＨ０は、最適化問題のコスト関数を表現したものであり、ハミルトニアンＨ０の基底状態が最適解となるようにＨ０を選択する。σはスピン変数であり、±１の値を取り得る。初期状態（時刻ｔ＝０）では、係数Γを非常に大きな値とし、時間の経過とともに係数Γを小さくし、最終的には０にする。最初は大きな量子揺らぎによって数多くの状態の重ね合わせを実現して状態探査をする。各時刻における瞬間的な基底状態を連続的にたどり、次第にΓが小さくなると、初期ハミルトニアンに比べてハミルトニアンＨ０の相対的な重みが大きくなり、最終的には、ハミルトニアンＨ０の基底状態に到達する。この状態で、最適化問題の解が得られ、所要の物理量を算出することができる。
　第１確率分布関数演算部１４は、横磁化ｍｘとスピンの横方向成分σｘの和の平均との差を変数とするデルタ関数と、磁場関数との乗算項を含む指数関数型演算子を用いて初期ハミルトニアンに対する確率分布関数を演算する。具体的には、第１確率分布関数演算部１４によって演算される確率分布関数は、式（５）のように表すことができる。
Figure JPOXMLDOC01-appb-I000004
　式（５）は、鈴木トロッター分解したときの確率分布を示す。式（５）において、βは絶対温度の逆数に比例する係数であり、τはトロッター数を示す。式（５）に示すように、磁場関数ｆ（ｍｘ）に関する効果を書き直すことによって、ｆ（ｍｘ）を含む、任意の量子揺らぎ（式（３）の右辺第２項で表した初期ハミルトニアンの問題）を、横磁場を持つ単純なものに変更することができる。
　式（５）に示すように、横磁化ｍｘとスピンのｘ方向成分σｘの和の平均との差を変数とするデルタ関数を導入することにより、横磁化ｍｘが、スピンのｘ方向成分σｘの和の平均と等しい場合には、デルタ関数が１となり、デルタ関数と磁場関数ｆ（ｍｘ）との乗算項は、結果として磁場関数ｆ（ｍｘ）に置き換えることができ、確率分布関数の指数関数型演算子からスピンのｘ方向成分（σｘ）の和の二次以上の高次項（ＸＸ相互作用を含む）を除くことができるので、負符号問題がなくなる。
　式（６）は、鈴木トロッター分解公式を示す。式（６）において、Ａ、Ｂは演算子であり、Ｌはトロッター数である。なお、本明細書では、トロッター数をτでも表している。
Figure JPOXMLDOC01-appb-I000005
　一般的に、量子系のハミルトニアンは、局所的な構成要素間の相互作用を表す局所ハミルトニアンの和で定義することができる。局所ハミルトニアン同士はお互いに非可換であるため、量子系のハミルトニアンの行列表現のサイズが大きくなり、計算コストが高くなる。そこで、式（６）で示す鈴木トロッター分解公式を用いて、指数演算子を行列表現のサイズが小さな局所ハミルトニアンの指数演算子の積に分解することができる。
　第２確率分布関数演算部１５は、第１確率分布関数演算部１４で演算した確率分布関数に含まれるデルタ関数を積分表示にして、磁場関数の導関数を含む指数関数型演算子を用いて系のハミルトニアンに対する確率分布関数を演算する。具体的には、第２確率分布関数演算部１５によって演算される確率分布関数は、式（７）のように表すことができる。
Figure JPOXMLDOC01-appb-I000006
　式（７）において、τはトロッター数であり、ｍチルダｘは、式（８）のように、磁場関数ｆ（ｍｘ）のｍｘについての導関数ｆ′（ｍｘ）で表すことができる。式（７）は、式（５）のデルタ関数を積分表示（フーリエ積分表示）することで、ｍチルダｘが導入されている。式（７）では、スピンのｘ方向成分σｘを持つハミルトニアンの問題を横磁場の問題として扱うことができ、鈴木トロッター分解を実行することによって、スピンのｘ方向成分σｘに関する項を式（９）のＢを係数として持つトロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。これにより、通常のコンピュータを用いてシミュレーションを実行することができる。なお、係数Ｂは、式（１０）で表すことができる。また、本明細書では、式（９）も第２確率分布関数演算部１５によって演算される確率分布関数と称する。式（９）において、Ｚ（ｍｘ）は規格化のための係数である。
　図２は鈴木トロッター分解の一例を示す模式図である。図２において、横軸は１次元格子上でのスピン変数を表し、いわゆる実空間方向を示す。縦軸は鈴木トロッター分解によって導入された方向（トロッター方向）であり、２次元の格子点上に状態変数が配置される。このように、鈴木トロッター分解によって量子モデルは、次元が一つ増えた状態空間を持つ古典モデルに変換されたと考えることができる。
　式（７）の指数関数部分に注目して、ｍｘについて極大値を求める。ｍｘについての極大を取るところだけに注目する鞍点法を利用することで、式（９）を導出することができる。その鞍点の条件（極大点となる）を表す式が、式（８）となり、式（９）でｍチルダｘが消えている。式（９）は、横磁化ｍｘが決まるとスピン変数σがどのような振る舞いをするかを示す。
　スピン変数更新部１６は、第１確率分布関数演算部１４によって演算して得られた確率分布に基づいて、複数のスピンそれぞれのスピン変数を更新する。より具体的は、スピン変数更新部１６は、第２確率分布関数演算部１５で演算して得られた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。スピン変数の更新は、例えば、スピン変数を選び、熱浴法又はメトロポリス法などの詳細釣り合いを満たすもので更新することができる。
　図３はスピン変数の更新の様子の一例を示す模式図である。左側の図は、現在のスピン変数の状態を示し、右図は任意のスピン変数として、格子ｉのスピン変数を選んだとする。なお、便宜上、格子ｉのスピン変数の周囲に４つのスピン変数があるとし、格子ｉのスピン変数しか変化できないとする（１スピンフリップ）。
　まず、左図の状態でのスピン変数を式（９）に代入して、現在の確率分布Ｐｐを算出する。次に、右図のように、選択した格子ｉのスピン変数を変更した状態でのスピン変数を式（９）に代入して、次の新しい確率分布Ｐｎを算出する。フリップ確率ＰｆをＰｆ＝Ｐｎ／Ｐｐにより算出する。一様乱数ｒを生成し、フリップ確率Ｐｆが乱数ｒより大きい場合、格子ｉのスピン変数を右図のとおり（スピンをフリップさせる）とし、フリップ確率Ｐｆが乱数ｒより大きくない場合、格子ｉのスピン変数を左図のとおり（スピンをフリップさせない）とする。上述の処理をすべての格子（すなわち、スピンの数Ｎ）について繰り返す。
　平衡状態判定部１７は、判定部としての機能を有し、更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。系が平衡状態になったか否かの判定は、例えば、系のエネルギーＥ、磁化ｍ（磁化の２乗、４乗などのｍ－メント）を算出（測定）して、変動がなければ系が平衡状態になったと判定することができる。系のエネルギーＥは、式（１）で算出する。磁化ｍは、全スピンについて和を計算し、その和をスピンの数で除算して求めることができる。
　第１磁化算出部１９は、平衡状態判定部１７で系が平衡状態になったと判定した場合、当該平衡状態での横磁化ｍｘを算出する。平衡状態での横磁化ｍｘは、式（１１）により算出した量の時間平均値で求めることができる。式（１１）において、ｉはスピンの場所（格子点）を示し、ｔはトロッター数を示し、τはトロッター総数を示し、Ｎはスピン総数を示す。
Figure JPOXMLDOC01-appb-I000007
　磁場算出部２０は、第１磁化算出部１９で算出した横磁化ｍｘ及び磁場関数ｆ（ｍｘ）に基づいて、複数のスピンに対する横磁場を算出する。横磁場（ｍチルダｘ）は、磁場関数ｆ（ｍｘ）を横磁化ｍｘで微分した式（８）に横磁化ｍｘを代入して算出することができる。
　磁場判定部２１は、磁場算出部２０で算出した横磁場（ｍチルダｘ）が定常状態であるか否かを判定する。横磁場が定常状態であるか否かは、平衡状態での横磁化ｍｘの値に応じて横磁場の値を適応的に変化させ、横磁場が変化しない場合、定常状態であると判定することができる。
　物理量算出部２２は、磁場判定部２１で横磁場が定常状態であると判定した場合、系に関する物理量を算出する。系に関する物理量は、例えば、エネルギーＥ、磁化ｍなどである。シミュレーションを続けながら、スピン変数に基づいて、エネルギーＥ、磁化ｍなどの物理量を繰り返し計算し、ある程度の時間が経過したときに物理量の時間平均を算出し、物理量の結果とする。時間平均によって任意の精度で計算することができ、時間を長くすれば、それに応じて精度を高めることができる。
　記憶部２４は、シミュレーションの必要なデータ、入力データ、シミュレーション中に得られた処理結果、出力データなどを記憶する。
　上述の構成により、一次相転移を回避しつつ負符号問題を解決して、シミュレーションを実施することができる。
　また、スピン変数更新部１６は、磁場判定部２１で横磁場が定常状態でないと判定した場合、系のハミルトニアンに対する確率分布関数に含まれる磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）を第１磁化算出部１９で算出した横磁化に基づいて更新し、更新した系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
　具体的には、式（９）の中にある係数Ｂが、式（１０）で示すように、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）を含んでいる。すなわち、ｍｘによってｆ′（ｍｘ）を更新すると、係数Ｂが変わり、結果として式（９）で演算される系の確率分布が変わる。
　図４は横磁化ｍｘと横磁場ｍチルダｘとの関係の一例を示す説明図である。ｍチルダｘは、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）である。例えば、磁場関数ｆ（ｍｘ）が式（４）で表される場合、導関数ｆ′（ｍｘ）は、横磁化ｍｘの一次関数となる。すなわち、横磁化ｍｘを変えることによって横磁場（ｍチルダｘ）を変えることができる。なお、従来の量子アニーリングでは、横磁場は一定であった。
　上述のように、横磁場が定常状態でない場合、解ではない別の横磁場又は横磁化が得られたことになる。そこで、系のハミルトニアンに対する確率分布を横磁化に応じて計算し、再度スピン変数を更新する処理を行うことにより、解に到達することができる。
　次に、本実施の形態のシミュレーション装置１００の動作について説明する。図５は本実施の形態のシミュレーション装置１００による適応的量子モンテカルロ法の処理手順の一例を示すフローチャートである。適応的量子モンテカルロ法での「適応的」とは、従来の量子アニーリング法を実現するための確率的手法である量子モンテカルロ法と区別するための用語であり、量子相移転（一次相転移）を回避しつつ負符号問題も回避することができ、通常のコンピュータにおいて実行可能なシミュレーション法であることを表す。なお、以下では、便宜上、処理の主体を制御部１０として説明する。
　制御部１０は、トロッター数、スピンの数を設定する（Ｓ１１）。トロッター数は、シミュレータ（コンピュータ）の性能にも依存するが、例えば、１２８とすることができる。スピンの数Ｎは、任意のサイズとすることができる。トロッター数もスピンの数も大きくすれば、計算精度を高めることができる。
　制御部１０は、スピン変数のシグマに初期値を設定する（Ｓ１２）。スピン変数のシグマに初期値を設定すると、横磁化ｍｘを計算することができるので、横磁場の値（ｍチルダｘ、すなわちｆ′（ｍｘ））を初期値に設定することができる。
　制御部１０は、スピン変数を乱数によって初期化し（Ｓ１３）、スピン変数を選び、熱浴法又はメトロポリス法（詳細釣り合いを満たすもの）で更新する（Ｓ１４）。これにより、系のハミルトニアンが目的ハミルトニアンＨ０の基底状態に向かって収束するようにする。
　制御部１０は、系が平衡状態になったか否かを判定する（Ｓ１５）。平衡状態になったか否かの判定は、エネルギーＥの値、磁化ｍなどが変動するか否かによって行うことができる。
　系が平衡状態でない場合（Ｓ１５でＮＯ）、制御部１０は、ステップＳ１４以降の処理を続ける。系が平衡状態である場合（Ｓ１５でＹＥＳ）、制御部１０は、平衡状態での横磁化ｍｘの値に応じて横磁場（ｍチルダｘ、すなわちｆ′（ｍｘ））の値を変更する（Ｓ１６）。
　制御部１０は、横磁場が変化するか、すなわち横磁場が定常状態であるか否かを判定する（Ｓ１７）。横磁場が変化する場合（Ｓ１７でＹＥＳ）、制御部１０は、解ではない別の横磁場又は横磁化が得られたとして、ステップＳ１４以降の処理を続ける。
　横磁場が変化しない場合（Ｓ１７でＮＯ）、制御部１０は、解が得られたとして物理量の測定（計算）を開始し（Ｓ１８）、所要の時間が経過したときに、測定量の時間平均を求め、物理量の結果とし、処理を終了する。
　図６は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。図６において、横軸は横磁場Γ（Gamma）を示し、縦軸はエネルギーＥを示す。図中Ｎはシミュレーション時のスピンの数を示し、exactγ＝１で示す実線は、手計算により求めた真の解（正解）を示す。図６から分かるように、シミュレーション結果は、真の解とほぼ一致し、特にスピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
　図７は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図であり、図８は本実施の形態の適応的量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。図７において、縦軸は磁化ｍを示す。また、図８において、縦軸は横磁化ｍｘを示す。図７及び図８においても、シミュレーション結果は、スピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
　次に、本実施の形態のシミュレーション装置１００が実行するテータ解析による量子モンテカルロ法について説明する。テータ解析による量子モンテカルロ法は、通常の量子モンテカルロ法である。
　制御部１０は、設定部としての機能を有し、予め横磁場の値を複数設定する。すなわち、横磁場の値を複数用意しておく。
　第２磁化算出部２３は、横磁場の値及び磁場関数の導関数の逆関数に基づいて横磁化を算出する。例えば、横磁場ｍチルダｘ＝ｆ′（ｍｘ）とすると、横磁化ｍｘは、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数で算出することができる。複数の横磁場に対して横磁化を算出することにより、横磁場と横磁化との関係をプロットすることができる。
　スピン変数更新部１６は、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数に制御部１０が設定した横磁場の値を割り当てた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
　平衡状態判定部１７は、スピン変数更新部１６で更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。
　第１磁化算出部１９は、系が平衡状態になった場合、当該平衡状態での横磁化を算出する。すなわち、予め容易した横磁場を用いて量子モンテカルロ法を実行して、横磁化を算出する。これにより、横磁場と横磁化との関係をプロットすることができる。
　物理量算出部２２は、第１磁化算出部１９及び第２磁化算出部２３で算出した横磁化が等しい場合、系に関する物理量を算出する。すなわち、磁場関数の導関数の逆関数に基づいて算出した横磁化と、量子モンテカルロ法を実行して計算された横磁化とが等しくなる場合、当該横磁化及び対応する横磁場を求めることにより、データ解析的なアプローチによっても解を求めることができる。
　図９は本実施の形態のシミュレーション装置１００が行うデータ解析による量子モンテカルロ法の処理手順の一例を示すフローチャートである。制御部１０は、横磁場（ｍチルダｘ、すなわち、ｆ′（ｍｘ））の値を複数設定し（Ｓ３１）、設定したそれぞれの横磁場の値を用いて量子モンテカルロ法を実行する（Ｓ３２）。
　制御部１０は、量子モンテカルロ法を実行することによって得られた横磁化の値と、当該横磁化が得られたときの横磁場とを対応付けて横磁場の値と横磁化の値との関係をプロットする（Ｓ３３）。なお、ここで、プロットするとは、実際にチャートに描く必要はなく、対応付けが分かるような形態であればよい。
　制御部１０は、磁場関数ｆ（ｍｘ）の導関数ｆ′（ｍｘ）の逆関数に基づいて横磁化を算出する（Ｓ３４）。具体的には、設定した横磁場の値それぞれを逆関数に代入して横磁化を算出することができる。
　制御部１０は、逆関数によって算出した横磁化が、プロット上の横磁化と一致するときの横磁化の値及び横磁場の値を特定する（Ｓ３５）。これにより、真の解が得られ、適応的量子モンテカルロ法と同様に物理量の結果を得ることができる。制御部１０は、処理を終了する。
　図１０は本実施のデータ解析による量子モンテカルロ法の概念を示す説明図である。図１０において、横軸は横磁場（ｍチルダｘ）の値を示し、縦軸は横磁化ｍｘの値を示す。図１０中、符号Ｐ１で示す曲線は、量子モンテカルロ法によって得られた結果をプロットしたものである。また、横磁場ｍチルダｘ＝ｆ′（ｍｘ）とすると、横磁化ｍｘは、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数で算出することができる。符号Ｐ２で示す直線は、ｍｘ＝ｆ′（ｍチルダｘ）の逆関数を満たす。そして、曲線Ｐ１と直線Ｐ２とが交差する点における横磁化ｍｘの値と横磁場ｍチルダｘの値が解となる。
　図１１は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第１例を示す説明図である。図１１において、横軸は横磁場Γ（Gamma）を示し、縦軸はエネルギーＥを示す。図中Ｎはシミュレーション時のスピンの数を示し、各Ｎに対応するチャートは、曲線Ｐ１と直線Ｐ２とが交差する点の集合である。exactγ＝１で示す実線は、手計算により求めた真の解（正解）を示す。図１１から分かるように、シミュレーション結果は、真の解とほぼ一致し、特にスピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
　図１２は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第２例を示す説明図であり、図１３は本実施の形態のデータ解析による量子モンテカルロ法によるシミュレーション結果の第３例を示す説明図である。図１２において、縦軸は磁化ｍを示す。また、図１３において、縦軸は横磁化ｍｘを示す。図１２及び図１３においても、シミュレーション結果は、スピン数Ｎを大きくすることにより、より真の解に近づくことが分かる。
　本実施の形態のシミュレーション装置１００は、ＣＰＵ（プロセッサ）、ＲＡＭ（メモリ）などを備えたコンピュータを用いて実現することもできる。すなわち、図５及び図９に示すような、各処理の手順を定めたコンピュータプログラムをコンピュータに備えられたＲＡＭ（メモリ）にロードし、コンピュータプログラムをＣＰＵ（プロセッサ）で実行することにより、シミュレーション装置１００を実現することができる。
　図１４は本実施の形態のシミュレーション装置の構成の他の例を示す説明図である。図１４において、符号３００は、通常のコンピュータである。コンピュータ３００は、制御部３０、入力部４０、出力部５０、外部Ｉ／Ｆ（インタフェース）部６０などを備える。制御部３０は、ＣＰＵ３１、ＲＯＭ３２、ＲＡＭ３３、Ｉ／Ｆ（インタフェース）３４などを備える。
　入力部４０は、シミュレーションのための入力データを取得する。出力部５０は、シミュレーション結果である出力データを出力する。Ｉ／Ｆ３４は、制御部３０と、入力部４０、出力部５０及び外部Ｉ／Ｆ部６０それぞれとの間のインタフェース機能を有する。
　外部Ｉ／Ｆ部６０は、コンピュータプログラムを記録した記録媒体Ｍ（例えば、ＤＶＤなどのメディア）からコンピュータプログラムを読み取ることが可能である。
　なお、図示していないが、記録媒体Ｍに記録されたコンピュータプログラムは、持ち運びが自由なメディアに記録されたものに限定されるものではなく、インターネット又は他の通信回線を通じて伝送されるコンピュータプログラムも含めることができる。また、コンピュータには、複数のプロセッサを搭載した１台のコンピュータ、あるいは、通信ネットワークを介して接続された複数台のコンピュータで構成されるコンピュータシステムも含まれる。
　上述のように、本実施の形態によれば、シミュレーションの計算時間が最適解を得ることができないほど長くなる量子相転移の問題を回避することができる。また、本実施の形態によれば、確率密度が負になってしまうという負符号問題も回避することができるので、限定的なモデルに限られることなく、幅広い範囲のモデルに対して、通常のコンピュータ上でシミュレーションを実行することが可能となり、量子モンテカルロ法の適用範囲を広げることができ、物質探索、あるいは設計のシミュレーションの範囲を拡大させることができる。また、本実施の形態のシミュレーション方法は、量子コンピュータの製作現場、人工知能、機械学習などの大規模な計算を必要とする技術分野で利用することができる。
　本実施の形態に係るシミュレーション装置は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部を備えることを特徴とする。
　本実施の形態に係るコンピュータプログラムは、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させることを特徴とする。
　本実施の形態に係るコンピュータでの読み取り可能な記録媒体は、コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、コンピュータに、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを実行させるためのコンピュータプログラムを記録したことを特徴とする。
　本実施の形態に係るシミュレーション方法は、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、算出した磁場が定常状態であるか否かを判定する処理と、前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理とを含むことを特徴とする。
　磁化演算部は、複数のスピンの所定方向成分の和の平均を当該所定方向の磁化として演算する。ｚ方向に向いたスピンを反転させる量子力学的な揺らぎとして、ｘ方向の横磁場を利用する場合、所定方向は、横方向であるｘ方向とすることができ、所定方向の磁化は、横方向の磁化ｍｘとすることができる。すなわち、横磁化ｍｘを、スピンのｘ方向成分のシグマの平均として演算する。
　初期ハミルトニアン演算部は、演算した磁化の一次項及び二次以上の項を含む磁場関数を初期ハミルトニアンとして演算する。磁場関数をｆ（ｍｘ）で表すと、磁場関数ｆ（ｍｘ）は、例えば、ｆ（ｍｘ）＝Γ・ｍｘ＋（γ／２）・（ｍｘ）2 と表すことができる。初期ハミルトニアンは、磁場関数で表すことができる。磁場関数に横磁化ｍｘの２乗の項を含めることにより、一次相転移の問題を回避することが可能となる。
　第１確率分布関数演算部は、演算した磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数と、磁場関数との乗算項を含む指数関数型演算子を用いて初期ハミルトニアンに対する確率分布関数を演算する。磁化とスピンの所定方向成分の和の平均との差を変数とするデルタ関数を導入することにより、横磁化ｍｘが、スピンの所定方向成分の和の平均と等しい場合には、デルタ関数が１となり、デルタ関数と磁場関数との乗算項は、結果として磁場関数に置き換えられるので、確率分布関数の指数関数型演算子からスピンの所定方向成分（σｘ）の和の２次以上の高次項を除くことができるので、負符号問題がなくなる。
　スピン変数更新部は、演算して得られた確率分布に基づいて、複数のスピンそれぞれのスピン変数を更新する。例えば、スピン変数を選び、熱浴法又はメトロポリス法などの詳細釣り合いを満たすもので更新することができる。
　判定部は、更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。系が平衡状態になったか否かの判定は、例えば、系のエネルギー、磁化（磁化の２乗、４乗などのｍ－メント）を算出（測定）して、変動がなければ系が平衡状態になったと判定することができる。
　第１磁化算出部は、系が平衡状態になった場合、当該平衡状態での所定方向の磁化（横磁化ｍｘ）を算出する。
　磁場算出部は、算出した磁化（横磁化）及び磁場関数に基づいて、複数のスピンに対する所定方向の磁場（横磁場）を算出する。横磁場は、磁場関数を横磁化で微分した式に横磁化を代入して算出することができる。
　磁場判定部は、算出した磁場（横磁場）が定常状態であるか否かを判定する。横磁場が定常状態であるか否かは、平衡状態での横磁化の値に応じて横磁場の値を適応的に変化させ、横磁場が変化しない場合、定常状態であると判定することができる。
　物理量算出部は、横磁場が定常状態であると判定した場合、系に関する物理量を算出する。系に関する物理量は、例えば、エネルギー、磁化などである。シミュレーションを続けながら、スピン変数に基づいて、エネルギー、磁化などの物理量を繰り返し計算し、ある程度の時間が経過したときに物理量の時間平均を算出し、物理量の結果とする。
　上述の構成により、一次相転移を回避しつつ負符号問題を解決して、シミュレーションを実施することができる。
　本実施の形態に係るシミュレーション装置は、前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、前記スピン変数更新部は、前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
　第２確率分布関数演算部は、第１確率分布関数演算部で演算した確率分布関数に含まれるデルタ関数を積分表示にして、磁場関数の導関数を含む指数関数型演算子を用いて系のハミルトニアンに対する確率分布関数を演算する。デルタ関数を積分表示にすることにより、磁場関数の導関数を導入することができ、鈴木トロッター分解を実行することで、スピンのｘ方向成分（σｘ）に関する項をトロッター間相互作用に置き換えることができ、数値計算を実行することが可能となる。
　スピン変数更新部は、第２確率分布関数演算部で演算して得られた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。これにより、通常のコンピュータを用いてシミュレーションを実行することができる。
　本実施の形態に係るシミュレーション装置は、前記スピン変数更新部は、前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする。
　スピン変数更新部は、横磁場が定常状態でないと判定した場合、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数を第１磁化算出部で算出した横磁化に基づいて更新した系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。
　横磁場が定常状態でない場合、解ではない別の横磁場又は横磁化が得られたことになる。そこで、系のハミルトニアンに対する確率分布を横磁化に応じて計算し、再度スピン変数を更新する処理を行うことにより、解に到達することができる。
　本実施の形態に係るシミュレーション装置は、予め前記所定方向の磁場の値を複数設定する設定部と、該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部とを備え、前記スピン変数更新部は、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、前記判定部は、前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、前記第１磁化算出部は、前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、前記物理量算出部は、前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする。
　設定部は、予め所定方向の磁場（横磁場）の値を複数設定する。すなわち、横磁場の値を複数用意しておく。
　第２磁化算出部は、横磁場の値及び磁場関数の導関数の逆関数に基づいて所定方向の磁化（横磁化）を算出する。複数の横磁場に対して横磁化を算出することにより、横磁場と横磁化との関係をプロットすることができる。
　スピン変数更新部は、系のハミルトニアンに対する確率分布関数に含まれる磁場関数の導関数に設定部で設定した磁場の値を割り当てた系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新する。判定部は、スピン変数更新部で更新したスピン変数に基づいて系が平衡状態になったか否かを判定する。第１磁化算出部は、系が平衡状態になった場合、当該平衡状態での横磁化を算出する。すなわち、予め容易した横磁場を用いて量子モンテカルロ法を実行して、横磁化を算出する。これにより、横磁場と横磁化との関係をプロットすることができる。
　物理量算出部は、第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、系に関する物理量を算出する。すなわち、磁場関数の導関数の逆関数に基づいて算出した横磁化と、量子モンテカルロ法を実行して計算された横磁化とが等しくなる場合、当該横磁化及び対応する横磁場を求めることにより、データ解析的なアプローチによっても解を求めることができる。
　１０、３０　制御部
　１１、４０　入力部
　１２　磁化演算部
　１３　初期ハミルトニアン演算部
　１４　第１確率分布関数演算部
　１５　第２確率分布関数演算部
　１６　スピン変数更新部
　１７　平衡状態判定部
　１８、５０　出力部
　１９　第１磁化算出部
　２０　磁場算出部
　２１　磁場判定部
　２２　物理量算出部
　２３　第２磁化算出部
　２４　記憶部
　３１　ＣＰＵ
　３２　ＲＯＭ
　３３　ＲＡＭ
　３４　Ｉ／Ｆ
　６０　外部Ｉ／Ｆ部

Claims (6)

  1. 　二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートするシミュレーション装置であって、
    　前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する磁化演算部と、
    　該磁化演算部で演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する初期ハミルトニアン演算部と、
    　前記磁化演算部で演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する第１確率分布関数演算部と、
    　該第１確率分布関数演算部で演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新するスピン変数更新部と、
    　該スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する判定部と、
    　該判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する第１磁化算出部と、
    　該第１磁化算出部で算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する磁場算出部と、
    　該磁場算出部で算出した磁場が定常状態であるか否かを判定する磁場判定部と、
    　該磁場判定部で前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する物理量算出部と
    　を備えることを特徴とするシミュレーション装置。
  2. 　前記第１確率分布関数演算部で演算した確率分布関数に含まれる前記デルタ関数を積分表示にして、前記磁場関数の導関数を含む指数関数型演算子を用いて前記系のハミルトニアンに対する確率分布関数を演算する第２確率分布関数演算部を備え、
    　前記スピン変数更新部は、
    　前記第２確率分布関数演算部で演算して得られた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする請求項１に記載のシミュレーション装置。
  3. 　前記スピン変数更新部は、
    　前記磁場判定部で前記磁場が定常状態でないと判定した場合、前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数を前記第１磁化算出部で算出した磁化に基づいて更新した前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新することを特徴とする請求項２に記載のシミュレーション装置。
  4. 　予め前記所定方向の磁場の値を複数設定する設定部と、
    　該設定部で設定した磁場の値及び前記磁場関数の導関数の逆関数に基づいて前記所定方向の磁化を算出する第２磁化算出部と
    　を備え、
    　前記スピン変数更新部は、
    　前記系のハミルトニアンに対する確率分布関数に含まれる前記磁場関数の導関数に前記設定部で設定した磁場の値を割り当てた前記系のハミルトニアンに対する確率分布に基づいて、スピン変数を更新し、
    　前記判定部は、
    　前記スピン変数更新部で更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定し、
    　前記第１磁化算出部は、
    　前記判定部で前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出し、
    　前記物理量算出部は、
    　前記第１磁化算出部及び第２磁化算出部で算出した磁化が等しい場合、前記系に関する物理量を算出することを特徴とする請求項３に記載のシミュレーション装置。
  5. 　コンピュータに、二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるコンピュータプログラムであって、
    　コンピュータに、
    　前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、
    　演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、
    　演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、
    　演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、
    　更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、
    　前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、
    　算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、
    　算出した磁場が定常状態であるか否かを判定する処理と、
    　前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理と
    　を実行させることを特徴とするコンピュータプログラム。
  6. 　二値を取り得る複数のスピンからなる系のハミルトニアンを初期ハミルトニアン及び目的ハミルトニアンで表し、初期状態において前記初期ハミルトニアンを大きな値に設定し、時間変化に伴って前記初期ハミルトニアンが前記目的ハミルトニアンに比べて小さくなるようにして前記系の平衡状態での物理量をシミュレートさせるシミュレーション方法であって、
    　前記複数のスピンの所定方向成分の和の平均を該所定方向の磁化として演算する処理と、
    　演算した磁化の一次項及び二次以上の項を含む磁場関数を前記初期ハミルトニアンとして演算する処理と、
    　演算した磁化と前記スピンの所定方向成分の和の平均との差を変数とするデルタ関数と、前記磁場関数との乗算項を含む指数関数型演算子を用いて前記初期ハミルトニアンに対する確率分布関数を演算する処理と、
    　演算して得られた確率分布に基づいて、前記複数のスピンそれぞれのスピン変数を更新する処理と、
    　更新したスピン変数に基づいて前記系が平衡状態になったか否かを判定する処理と、
    　前記系が平衡状態になったと判定した場合、前記平衡状態での前記所定方向の磁化を算出する処理と、
    　算出した磁化及び前記磁場関数に基づいて、前記複数のスピンに対する前記所定方向の磁場を算出する処理と、
    　算出した磁場が定常状態であるか否かを判定する処理と、
    　前記磁場が定常状態であると判定した場合、前記系に関する物理量を算出する処理と
    　を含むことを特徴とするシミュレーション方法。
     
PCT/JP2017/023378 2016-10-20 2017-06-26 シミュレーション装置、コンピュータプログラム及びシミュレーション方法 WO2018074006A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CA3041221A CA3041221A1 (en) 2016-10-20 2017-06-26 Simulation device, computer program, and simulation method
US16/341,024 US20190235033A1 (en) 2016-10-20 2017-06-26 Simulation Device, Recording Medium, and Simulation Method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2016206265A JP2018067200A (ja) 2016-10-20 2016-10-20 シミュレーション装置、コンピュータプログラム及びシミュレーション方法
JP2016-206265 2016-10-20

Publications (1)

Publication Number Publication Date
WO2018074006A1 true WO2018074006A1 (ja) 2018-04-26

Family

ID=62019125

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/023378 WO2018074006A1 (ja) 2016-10-20 2017-06-26 シミュレーション装置、コンピュータプログラム及びシミュレーション方法

Country Status (4)

Country Link
US (1) US20190235033A1 (ja)
JP (1) JP2018067200A (ja)
CA (1) CA3041221A1 (ja)
WO (1) WO2018074006A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020100698A1 (ja) * 2018-11-12 2020-05-22 国立大学法人京都大学 シミュレーション装置、コンピュータプログラム及びシミュレーション方法
WO2020245877A1 (ja) * 2019-06-03 2020-12-10 日本電気株式会社 量子アニーリング計算装置、量子アニーリング計算方法および量子アニーリング計算プログラム
WO2021085297A1 (ja) * 2019-10-29 2021-05-06 国立大学法人東北大学 組合せ最適化問題処理装置、組合せ最適化問題処理方法及びプログラム

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109190247B (zh) * 2018-09-01 2020-11-06 刘照森 优化的量子蒙特-卡洛模拟方法在研究复杂磁系统中的应用
JP7228425B2 (ja) * 2019-03-19 2023-02-24 株式会社東芝 計算装置、表示装置およびプログラム
JP7171520B2 (ja) 2019-07-09 2022-11-15 株式会社日立製作所 機械学習システム
JP7341804B2 (ja) * 2019-09-06 2023-09-11 株式会社日立製作所 情報処理装置および情報処理方法
JP2021144622A (ja) 2020-03-13 2021-09-24 富士通株式会社 最適化装置、最適化プログラム、および最適化方法
EP3982301A1 (en) 2020-10-09 2022-04-13 Fujitsu Limited Apparatus, program, and method for optimization
CN114528996B (zh) * 2022-01-27 2023-08-04 本源量子计算科技（合肥）股份有限公司 一种目标体系试验态初始参数的确定方法、装置及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015103375A1 (en) * 2014-01-06 2015-07-09 Google Inc. Constructing and programming quantum hardware for quantum annealing processes
US20160071021A1 (en) * 2014-09-09 2016-03-10 D-Wave Systems Inc. Systems and methods for improving the performance of a quantum processor via reduced readouts

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015103375A1 (en) * 2014-01-06 2015-07-09 Google Inc. Constructing and programming quantum hardware for quantum annealing processes
US20160071021A1 (en) * 2014-09-09 2016-03-10 D-Wave Systems Inc. Systems and methods for improving the performance of a quantum processor via reduced readouts

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
HIDETOSHI NISHIMORI: "Ryoshi Annealing no Suri", BUSSEI KENKYU·DENSHI-BAN, vol. 3, no. 3, August 2014 (2014-08-01), pages 1 - 24 *
SEI SUZUKI: "Yoko Sogo Sayo Annealing", GRANT- IN-AID FOR SCIENTIFIC RESEARCH ON PRIORITY AREAS 'WORKSHOP ON DEEPENING AND EXPANSION OF STATISTICAL MECHANICAL INFORMATICS, 18 December 2006 (2006-12-18), pages 1 - 10 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020100698A1 (ja) * 2018-11-12 2020-05-22 国立大学法人京都大学 シミュレーション装置、コンピュータプログラム及びシミュレーション方法
WO2020245877A1 (ja) * 2019-06-03 2020-12-10 日本電気株式会社 量子アニーリング計算装置、量子アニーリング計算方法および量子アニーリング計算プログラム
WO2021085297A1 (ja) * 2019-10-29 2021-05-06 国立大学法人東北大学 組合せ最適化問題処理装置、組合せ最適化問題処理方法及びプログラム

Also Published As

Publication number Publication date
JP2018067200A (ja) 2018-04-26
US20190235033A1 (en) 2019-08-01
CA3041221A1 (en) 2018-04-26

Similar Documents

Publication Publication Date Title
WO2018074006A1 (ja) シミュレーション装置、コンピュータプログラム及びシミュレーション方法
Pandey et al. Machine learning based surrogate modeling approach for mapping crystal deformation in three dimensions
Li et al. State estimation on positive Markovian jump systems with time-varying delay and uncertain transition probabilities
Liu et al. A highly efficient and accurate exponential semi-implicit scalar auxiliary variable (ESI-SAV) approach for dissipative system
Xu et al. Kernel-based random vector functional-link network for fast learning of spatiotemporal dynamic processes
Bhattacharyya et al. A LATIN-based model reduction approach for the simulation of cycling damage
CN107016155A (zh) 非线性pde和线性求解器的收敛估计
US11150615B2 (en) Optimization device and control method of optimization device
Yu et al. Complex dynamics in biological systems arising from multiple limit cycle bifurcation
Sarkar et al. Self-learning emulators and eigenvector continuation
Deb et al. An optimality theory based proximity measure for evolutionary multi-objective and many-objective optimization
WO2020100698A1 (ja) シミュレーション装置、コンピュータプログラム及びシミュレーション方法
Lötzsch et al. Learning the solution operator of boundary value problems using graph neural networks
Amalinadhi et al. On physics-informed deep learning for solving navier-stokes equations
Chen et al. Gpt-pinn: Generative pre-trained physics-informed neural networks toward non-intrusive meta-learning of parametric pdes
Saito et al. An iterative improvement method for HHL algorithm for solving linear system of equations
Yao et al. Multiadam: Parameter-wise scale-invariant optimizer for multiscale training of physics-informed neural networks
Danciu A CNN-based approach for a class of non-standard hyperbolic partial differential equations modeling distributed parameters (nonlinear) control systems
Ding et al. Functional order-reduced Gaussian Processes based machine-learning emulators for probabilistic constitutive modelling
JP6651254B2 (ja) シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置
Chandra et al. Neural oscillators for magnetic hysteresis modeling
Rothkopf et al. A symmetry and Noether charge preserving discretization of initial value problems
Sheu et al. Using the combination of GM (1, 1) and Taylor approximation method to predict the academic achievement of student
Dwight et al. Adaptive uncertainty quantification for computational fluid dynamics
Tian et al. M-matrix-based state observer design for genetic regulatory networks with mixed delays

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17861661

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 3041221

Country of ref document: CA

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17861661

Country of ref document: EP

Kind code of ref document: A1