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

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

Info

Publication number
JP2018067200A
JP2018067200A JP2016206265A JP2016206265A JP2018067200A JP 2018067200 A JP2018067200 A JP 2018067200A JP 2016206265 A JP2016206265 A JP 2016206265A JP 2016206265 A JP2016206265 A JP 2016206265A JP 2018067200 A JP2018067200 A JP 2018067200A
Authority
JP
Japan
Prior art keywords
magnetic field
magnetization
hamiltonian
function
probability distribution
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
JP2016206265A
Other languages
English (en)
Inventor
真之 大関
Masayuki Ozeki
真之 大関
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kyoto University
Original Assignee
Kyoto University
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 Kyoto University filed Critical Kyoto University
Priority to JP2016206265A priority Critical patent/JP2018067200A/ja
Priority to CA3041221A priority patent/CA3041221A1/en
Priority to US16/341,024 priority patent/US20190235033A1/en
Priority to PCT/JP2017/023378 priority patent/WO2018074006A1/ja
Publication of JP2018067200A publication Critical patent/JP2018067200A/ja
Ceased legal-status Critical Current

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

Abstract

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

Description

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

Claims (6)

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

Priority Applications (4)

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

Applications Claiming Priority (1)

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

Publications (1)

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

Family

ID=62019125

Family Applications (1)

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

Country Status (4)

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

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020189315A1 (ja) * 2019-03-19 2020-09-24 株式会社東芝 計算装置、表示装置およびプログラム
EP3879463A1 (en) 2020-03-13 2021-09-15 Fujitsu Limited Optimization apparatus, optimization program, and optimization method
EP3982301A1 (en) 2020-10-09 2022-04-13 Fujitsu Limited Apparatus, program, and method for optimization
US11715036B2 (en) 2019-07-09 2023-08-01 Hitachi, Ltd. Updating weight values in a machine learning system

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109190247B (zh) * 2018-09-01 2020-11-06 刘照森 优化的量子蒙特-卡洛模拟方法在研究复杂磁系统中的应用
JP2020080006A (ja) * 2018-11-12 2020-05-28 国立大学法人京都大学 シミュレーション装置、コンピュータプログラム及びシミュレーション方法
WO2020245877A1 (ja) * 2019-06-03 2020-12-10 日本電気株式会社 量子アニーリング計算装置、量子アニーリング計算方法および量子アニーリング計算プログラム
JP7341804B2 (ja) * 2019-09-06 2023-09-11 株式会社日立製作所 情報処理装置および情報処理方法
AU2020376359A1 (en) * 2019-10-29 2022-05-12 Tohoku University Combinatorial optimization problem processing device, combinatorial optimization problem processing method, and program
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
西森 秀稔: "量子アニーリングの数理", 物性研究・電子版, vol. 3, no. 3, JPN6017026175, August 2014 (2014-08-01), pages 1 - 24, ISSN: 0004383715 *
鈴木 正: "横相互作用量子アニーリング", 科研費特定領域研究「情報統計力学の深化と展開」平成18年度研究成果発表会, JPN6017026176, 18 December 2006 (2006-12-18), pages 1 - 10, ISSN: 0004383716 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020189315A1 (ja) * 2019-03-19 2020-09-24 株式会社東芝 計算装置、表示装置およびプログラム
JP2020154524A (ja) * 2019-03-19 2020-09-24 株式会社東芝 計算装置、表示装置およびプログラム
US11715036B2 (en) 2019-07-09 2023-08-01 Hitachi, Ltd. Updating weight values in a machine learning system
EP3879463A1 (en) 2020-03-13 2021-09-15 Fujitsu Limited Optimization apparatus, optimization program, and optimization method
US11714936B2 (en) 2020-03-13 2023-08-01 Fujitsu Limited Optimization apparatus, optimization method, and non-transitory computer-readable storage medium for storing optimization program
EP3982301A1 (en) 2020-10-09 2022-04-13 Fujitsu Limited Apparatus, program, and method for optimization

Also Published As

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

Similar Documents

Publication Publication Date Title
WO2018074006A1 (ja) シミュレーション装置、コンピュータプログラム及びシミュレーション方法
Umar et al. Intelligent computing for numerical treatment of nonlinear prey–predator models
Wang et al. Fractional-order gradient descent learning of BP neural networks with Caputo derivative
Dehghannasiri et al. Optimal experimental design for materials discovery
Kaufman et al. CALPHAD, first and second generation–Birth of the materials genome
Noel A new gradient based particle swarm optimization algorithm for accurate computation of global minimum
CORNUET et al. Adaptive multiple importance sampling
Yang et al. Parameter identification of the generalized Prandtl–Ishlinskii model for piezoelectric actuators using modified particle swarm optimization
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
US11150615B2 (en) Optimization device and control method of optimization device
CN107016155A (zh) 非线性pde和线性求解器的收敛估计
Bhattacharyya et al. A LATIN-based model reduction approach for the simulation of cycling damage
Waseem et al. Investigation of singular ordinary differential equations by a neuroevolutionary approach
Deb et al. An optimality theory based proximity measure for evolutionary multi-objective and many-objective optimization
WO2020100698A1 (ja) シミュレーション装置、コンピュータプログラム及びシミュレーション方法
Ayad et al. Parametric analysis for genetic algorithms handling parameters
Lötzsch et al. Learning the solution operator of boundary value problems using graph neural networks
Lu et al. Adaptive online data-driven closed-loop parameter control strategy for swarm intelligence algorithm
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
Santana et al. A local mesh free numerical method with automatic parameter optimization
Matsko et al. Adaptive fuzzy decision tree with dynamic structure for automatic process control system o of continuous-cast billet production
JP6651254B2 (ja) シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20161102

A80 Written request to apply exceptions to lack of novelty of invention

Free format text: JAPANESE INTERMEDIATE CODE: A80

Effective date: 20161102

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20191011

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20201110

A045 Written measure of dismissal of application [lapsed due to lack of payment]

Free format text: JAPANESE INTERMEDIATE CODE: A045

Effective date: 20210323