(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023165119
(43)【公開日】2023-11-15
(54)【発明の名称】解探索プログラム、解探索方法、および情報処理装置
(51)【国際特許分類】
G06N 10/60 20220101AFI20231108BHJP
【ＦＩ】
G06N10/60
【審査請求】未請求
【請求項の数】10
【出願形態】ＯＬ
(21)【出願番号】P 2022075774
(22)【出願日】2022-05-02
(71)【出願人】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】110002918
【氏名又は名称】弁理士法人扶桑国際特許事務所
(72)【発明者】
【氏名】櫛部 大介
(57)【要約】
【課題】量子アニーリングにおける計算量を削減する。
【解決手段】情報処理装置１０は、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化を求めることでイジングモデルの基底状態を求解する。その際、情報処理装置１０は、イジングモデルの初期状態から、第１のアニーリングスケジュールに従って磁場の強度を低減させる。そして情報処理装置１０は、イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って磁場の強度を低減させる。
【選択図】図１
【解決手段】情報処理装置１０は、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化を求めることでイジングモデルの基底状態を求解する。その際、情報処理装置１０は、イジングモデルの初期状態から、第１のアニーリングスケジュールに従って磁場の強度を低減させる。そして情報処理装置１０は、イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って磁場の強度を低減させる。
【選択図】図１
【特許請求の範囲】
【請求項1】
求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化を求めることで前記イジングモデルの基底状態を求解するとき、前記イジングモデルの初期状態から、第１のアニーリングスケジュールに従って前記磁場の強度を低減させ、前記イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って前記磁場の強度を低減させる、
処理をコンピュータに実行させる解探索プログラム。
処理をコンピュータに実行させる解探索プログラム。
【請求項2】
前記第２のアニーリングスケジュールは、前記イジングモデルにおいて量子相転移が発生した後に前記第１のアニーリングスケジュールに従って前記磁場の強度を低減させた場合よりも、前記磁場の最小値に至るまでの前記磁場の更新回数が少ない、
請求項１記載の解探索プログラム。
請求項１記載の解探索プログラム。
【請求項3】
所定のパラメータの値が異なる複数の求解条件それぞれについて前記イジングモデルの基底状態を求解し、求解した結果に基づいて、前記量子相転移が発生するときの前記イジングモデルに印加される前記磁場の強さが所定の範囲内となる前記パラメータの値を決定する、
処理を前記コンピュータにさらに実行させる請求項１または２に記載の解探索プログラム。
処理を前記コンピュータにさらに実行させる請求項１または２に記載の解探索プログラム。
【請求項4】
前記パラメータの値を決定する処理では、前記パラメータの値が異なる複数の求解条件それぞれについての、最も低いエネルギー状態のスピン配置が得られる頻度に基づいて、前記量子相転移が発生するときの前記イジングモデルに印加される前記磁場の強さが所定の範囲内となる前記パラメータの値を決定する、
請求項３記載の解探索プログラム。
請求項３記載の解探索プログラム。
【請求項5】
印加される前記磁場を低減させたときの前記イジングモデルの状態変化を求める過程における各状態でのエネルギーの変化に基づいて量子相転移の発生の有無を判断する、
処理を前記コンピュータにさらに実行させる請求項１記載の解探索プログラム。
処理を前記コンピュータにさらに実行させる請求項１記載の解探索プログラム。
【請求項6】
前記イジングモデルのハミルトニアンは、指定された次数の非線形の外部磁場を示す項を含んでいる、
請求項１記載の解探索プログラム。
請求項１記載の解探索プログラム。
【請求項7】
前記イジングモデルの基底状態を求解するとき、前記磁場の強度を低減させていく実数時間発展の第１処理と、前記第１処理で前記磁場の強度を更新するごとに、虚数時間発展法に基づき、前記イジングモデルのエネルギーを低減させる第２処理とを実行する、
請求項１または６記載の解探索プログラム。
請求項１または６記載の解探索プログラム。
【請求項8】
前記第２処理の実行後に、前記第１処理により前記磁場の強度を変更するとき、直前に実行した前記第２処理の終状態を、前記磁場の強度を更新後の前記第１処理の初期状態とし、前記終状態と前記初期状態とのエネルギー差が所定の閾値より小さくなる範囲で、前記磁場の強度を変更する、
請求項７記載の解探索プログラム。
請求項７記載の解探索プログラム。
【請求項9】
求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化を求めることで前記イジングモデルの基底状態を求解するとき、前記イジングモデルの初期状態から、第１のアニーリングスケジュールに従って前記磁場の強度を低減させ、前記イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って前記磁場の強度を低減させる、
処理をコンピュータに実行させる解探索方法。
処理をコンピュータに実行させる解探索方法。
【請求項10】
求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化を求めることで前記イジングモデルの基底状態を求解するとき、前記イジングモデルの初期状態から、第１のアニーリングスケジュールに従って前記磁場の強度を低減させ、前記イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って前記磁場の強度を低減させる処理部、
を有する情報処理装置。
を有する情報処理装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、解探索プログラム、解探索方法、および情報処理装置に関する。
【背景技術】
【0002】
自然科学や社会科学において頻出する問題として、評価関数（エネルギー関数とも呼ばれる）の最小値（または最小値を与える評価関数の状態変数の値の組合せ）を求める最小値求解問題（または組合せ最適化問題と呼ばれる）がある。近年、このような問題を磁性体のスピンの振る舞いを表すモデルであるイジングモデルで定式化する動きが加速している。
【0003】
よく使用される解法は、マルコフ連鎖モンテカルロ（ＭＣＭＣ：Markov Chain Monte Carlo）法を用い確率過程を導入し、ボルツマン分布など特定の分布を定義した上で温度を導入し、温度をシミュレーション的に下げる方法である。この解法はシミュレーテッドアニーリング（ＳＡ：Simulated Annealing）法と呼ばれる。しかし、ＳＡ法で厳密解に到達するためには温度スケジュールを時間に関して対数の逆数でとることとなる。すなわち、高速に温度の冷却を行うと、必ずしも解が最適解に達しない。
【0004】
このような状況下で、近年量子力学に基づく計算技術および計算機の開発の動きが活発化している。この動きの技術的な基盤は量子アニーリング型量子コンピュータの実現である。量子アニーリング型の量子コンピュータは量子アニーリング（ＱＡ：Quantum Annealing）法が理論的な基盤になっている。量子コンピュータは普段使用している電子回路をベースにした古典的なコンピュータでは現実的な計算時間内で解くことが難しい問題を解くことができると期待されている。現段階ではすべての問題が量子コンピュータで高速化できるわけではないが、一部の問題が量子コンピュータで高速に解けると期待されている。代表的な問題がＲＳＡ暗号など暗号通信の基礎になっている素因数分解の問題である。このように、量子アニーリング型ではあるが量子コンピュータが実現されたため、量子コンピュータの応用範囲を探索する動きが加速している。
【0005】
ＱＡにおいては、横磁場と呼ばれる磁場を加えた項が初期状態として用意される。なお外部磁場を加えることを掃引と呼ぶこともある。ここから励起状態に励起させないように磁場を弱めることで、求めたい最適化問題のハミルトニアンに十分ゆっくり遷移させることができる。このようにして求解対象の量子イジングハミルトニアンの基底状態を求めることができる。
【0006】
一方でＱＡには、系のハミルトニアンによっては基底状態と励起状態のエネルギーギャップが小さくなり、断熱遷移にかける時間が伸びてしまうという問題がある。具体的には断熱遷移をしていく過程で、中間のハミルトニアンにおける基底状態と第一励起状態のエネルギー差が小さくなってしまう問題である。この問題を解決するため、時間依存シュレーディンガー方程式を虚数時間で解くことにより、現象論的なエネルギー散逸効果を取り入れた計算法も提案されている。
【先行技術文献】
【特許文献】
【0007】
【特許文献1】特開２０２１－１４４６２２号公報
【特許文献2】特開２０２２－０６２７６０号公報
【非特許文献】
【0008】
【非特許文献1】Hidetoshi Nishimori and Kabuki Takada, “Exponential Enhancement of the Efficiency of Quantum Annealing by Non-Stoquastic Hamiltonians”, Frontiers in ICT, Volume 4, Article 2, February 2017
【非特許文献2】T. Albash and D. A. Lidar, “Demonstration of a scaling advantage for a quantum annealer over simulated annealing”, arXiv:1705.07452v3, 6 August 2018
【非特許文献3】S. Morita, and H. Nishimori, ”Mathematical foundation of quantum annealing”, J. Math. Phys. 49, 125210(2008), 15 December 2008
【発明の概要】
【発明が解決しようとする課題】
【0009】
量子アニーリングの過程で量子相転移が発生し、計算対象のイジングモデルが励起状態になる場合がある。励起状態になった後は、基底状態に落ちずに励起状態のまま収束する。そのため量子相転移の発生後は横磁場の強度をゆっくり弱めても、イジングモデルの基底状態を得ることはできない。しかし、従来は、量子相転移が発生した後も、発生前と同じアニーリングスケジュールに従って横磁場が弱められており、余分な処理が発生している。
【0010】
１つの側面では、本件は、量子アニーリングにおける計算量を削減することを目的とする。
【課題を解決するための手段】
【0011】
１つの案では、以下の処理をコンピュータに実行させる解探索プログラムが提供される。
コンピュータは、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化を求めることでイジングモデルの基底状態を求解するとき、イジングモデルの初期状態から、第１のアニーリングスケジュールに従って磁場の強度を低減させ、イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って磁場の強度を低減させる。
コンピュータは、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化を求めることでイジングモデルの基底状態を求解するとき、イジングモデルの初期状態から、第１のアニーリングスケジュールに従って磁場の強度を低減させ、イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って磁場の強度を低減させる。
【発明の効果】
【0012】
１態様によれば、量子アニーリングにおける計算量を削減することができる。
【図面の簡単な説明】
【0013】
【図1】第１の実施の形態に係る解探索方法の一例を示す図である。
【図2】第２の実施の形態のシステム構成の一例を示す図である。
【図3】古典コンピュータのハードウェアの一例を示す図である。
【図4】古典コンピュータが有する機能を示すブロック図である。
【図5】虚数時間発展法を用いた熱的散逸の模式図である。
【図6】４都市のＴＳＰにおける都市の配置例を示す図である。
【図7】初期状態にランダムに付加したノイズに応じた計算結果の違いを示す図である。
【図8】８都市ＴＳＰの一例を示す図である。
【図9】８都市ＴＳＰの解探索結果の一例を示す図である。
【図10】Ｇ2の関数として得られた系の全エネルギー値を示す図である。
【図11】外部磁場の時間変化の一例を示す図である。
【図12】解析条件設定画面の一例を示す図である。
【図13】スピンの磁気の強さ決定処理の手順の一例を示すフローチャートである。
【図14】ｔqptの決定を伴うサンプリング処理の手順の一例を示すフローチャートの前半である。
【図15】ｔqptの決定を伴うサンプリング処理の手順の一例を示すフローチャートの後半である。
【図16】（１．１）式型の横磁場項と（２．１）式型の非線形横磁場を含んだ場合のＧ2に関する断熱ポテンシャルを示す図である。
【図17】横磁場の強さＧ2に関するポテンシャルエネルギーの一例を示す図である。
【図18】断熱ポテンシャルが交差する部分を拡大した図である。
【図19】非線形項として２次の横磁場を入れたハミルトニアン（２．１）式を対角化した例を示す図である。
【図20】４都市ＴＳＰの場合のエネルギー準位と第ｎ励起状態のエネルギーＥnと基底状態のエネルギーＥ0との差を示す図である。
【図21】Ｇ2の弱め方の手順の一例を示すフローチャートである。
【図22】４都市ＴＳＰのハミルトニアンを数値的に対角化したときの固有ベクトルの一例を示す図である。
【図23】固有ベクトルのインデックスと固有状態との対応関係の一例を示す図である。
【図24】量子相転移の検出を伴うアニーリング処理の手順の一例を示すフローチャートである。
【発明を実施するための形態】
【0014】
以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
第１の実施の形態は、量子アニーリングによるイジングモデルの基底状態の探索における計算量を削減する解探索方法である。
〔第１の実施の形態〕
第１の実施の形態は、量子アニーリングによるイジングモデルの基底状態の探索における計算量を削減する解探索方法である。
【0015】
図１は、第１の実施の形態に係る解探索方法の一例を示す図である。図１には、解探索方法を実施する情報処理装置１０が示されている。情報処理装置１０は、例えば解探索プログラムを実行することにより、解探索方法を実施することができる。
【0016】
情報処理装置１０は、例えば量子アニーリング型量子コンピュータ１に接続されている。情報処理装置１０は、量子アニーリング型量子コンピュータ１を制御し、量子アニーリングを用いて、求解対象の問題に応じたイジングモデルの基底状態を探索することで、例えば組合せ最適化問題の解を計算する。
【0017】
情報処理装置１０は、記憶部１１と処理部１２とを有する。記憶部１１は、例えば情報処理装置１０が有するメモリまたはストレージ装置である。処理部１２は、例えば情報処理装置１０が有するプロセッサまたは演算回路である。
【0018】
記憶部１１には、例えば求解対象の問題を表すイジングモデルの情報が格納される。イジングモデルの情報には、そのイジングモデルに対応するハミルトニアンの式が含まれる。
【0019】
処理部１２は、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化を求めることでイジングモデルの基底状態を求解する。その際、処理部１２は、イジングモデルの初期状態から、第１のアニーリングスケジュールに従って磁場の強度を低減させる。そして処理部１２は、イジングモデルにおいて量子相転移が発生した後は、第２のアニーリングスケジュールに従って磁場の強度を低減させる。
【0020】
このように、量子相転移前後において、アニーリングのスケジュールを変更することで、量子相転移の発生前と発生後とのそれぞれに応じた適切なアニーリングスケジュールを適用することができる。その結果、計算量の削減が可能となる。
【0021】
すなわち量子相転移発生前と量子相転移発生後とでは、イジングモデルの量子状態の変わりやすさが異なるため、共通のアニーリングスケジュールを適用するのは不適切である。すなわち量子相転移発生後は、イジングモデルの量子状態が他の状態に遷移しない。そのため、量子相転移が発生した後は、磁場の低減速度を速めても、解の探索精度に影響がない。
【0022】
そこで処理部１２は、例えば、第２のアニーリングスケジュールでは、イジングモデルにおいて量子相転移が発生した後に第１のアニーリングスケジュールに従って磁場の強度を低減させた場合よりも、磁場の最小値に至るまでの磁場の更新回数を少なくする。これにより、量子相転移が発生した後は、急速に磁場の強度が低減されることとなり、磁場が更新されるごとに実行される量子試行の回数が削減される。その結果、計算量も削減され、解探索時間も短縮される。
【0023】
また処理部１２は、量子相転移が発生するときのイジングモデルに印加される磁場の強さが所定の範囲内の値となるように、イジングモデルに含まれる所定のパラメータの値を決定することもできる。例えば処理部１２は、所定のパラメータの値が異なる複数の求解条件それぞれについてイジングモデルの基底状態を求解する。そして処理部１２は、求解した結果に基づいて、量子相転移が発生するときのイジングモデルに印加される磁場の強さが所定の範囲内となるパラメータの値を決定する。例えば処理部１２は、パラメータの値が異なる複数の求解条件それぞれについての、最も低いエネルギーの状態のスピン配置が得られる頻度に基づいて、パラメータの値を決定する。例えば最も低いエネルギーの状態のスピン配置が得られる頻度が最も高くなるパラメータの値が、量子相転移が発生するときのイジングモデルに印加される磁場の強さが所定の範囲内となるパラメータの値として決定される。
【0024】
処理部１２は、決定したパラメータの値をイジングモデルに設定して、イジングモデルの基底状態の探索を行う。これにより、量子相転移が発生するときの磁場の強さを適切な範囲に制御することができる。例えば量子相転移が発生する磁場の強さが大きすぎると、解探索の精度が低下する。逆に、量子相転移が発生する磁場の強さが小さすぎると、計算時間が長期化する。量子相転移が発生する磁場の強さを適切な範囲に制御できることで、要求する解探索精度に応じた効率的な解探索が可能となる。
【0025】
なお、量子相転移が発生する磁場の強さに影響を与えるパラメータとしては、例えばスピンの磁気の強さがある。スピンの磁気の強さは、イジングモデルに印可する磁場を示す横磁場項に含まれるパラメータである。
【0026】
なお量子相転移が発生するときの磁場の強さ、または量子アニーリングにおける時刻が分かっていれば、処理部１２にその値を設定しておくことで、計画的にアニーリングスケジュールの変更が可能である。量子相転移が発生するときの磁場の強さ、または量子アニーリングにおける時刻が不明であれば、量子アニーリングの過程で、量子相転移の発生の有無を検知することも可能である。
【0027】
例えば処理部１２は、印加される磁場を低減させたときのイジングモデルの状態変化を求める過程の各状態で得られるエネルギーの変化に基づいて量子相転移の発生の有無を判断する。量子相転移の発生の有無は、例えばエネルギーを横磁場の強さなどによる２回微分を求め、その２回微分の符号の変化などから判断可能である。
【0028】
また処理部１２は、イジングモデルのハミルトニアンに、指定された次数の非線形の外部磁場を示す項を含めることもできる。非線形の外部磁場を加えることにより、励起状態への遷移を抑制することができる。例えば処理部１２は、非線形横磁場の強度を量子アニーリングスケジュールに応じて変更する。これにより、基底状態が得られる確率を向上させることができる。
【0029】
さらに処理部１２は、イジングモデルの基底状態の求解において、磁場の強度を低減させていく実数時間発展の第１処理と、虚数時間発展法に基づき、イジングモデルのエネルギーを低減させる第２処理とを実行することができる。第２の処理は、例えば第１処理で磁場の強度を更新するごとに実行される。虚数時間発展法を用いることで熱散逸を発生させることができる。この熱散逸機構の導入により、断熱遷移条件を厳密に守らずとも基底状態へ緩和できるようになる。
【0030】
なお虚数時間発展法による熱散逸を発生させても、イジングモデルの状態が励起状態に遷移してしまう場合がある。励起状態に遷移してしまうと、その後、虚数時間発展を行っても基底状態に到達できない。そこで処理部１２は、イジングモデルのハミルトニアンに非線形の外部磁場を加えた上で、虚数時間発展法による熱散逸を実施してもよい。これにより、基底状態に遷移する確率を向上させることができる。
【0031】
なお処理部１２は、イジングモデルの初期状態として、異なるノイズを付加した複数の初期状態から解探索を実施することで、系の波動関数に多様性を持たせてもよい。これにより、イジングモデルの基底状態を確実に求解することが可能となる。
【0032】
虚数時間発展法に係る第２処理を適用する場合、処理部１２は、外部磁場を更新する際、磁場の更新前後の状態間でのエネルギーの差が所定の閾値以下となるように、外部磁場の変化量を制御してもよい。例えば処理部１２は、虚数時間発展法に係る第２処理の実行後に、第１処理により磁場の強度を変更するとき、直前に実行した第２処理の終状態と、磁場の強度を更新後の第１処理の初期状態とのエネルギーを計算する。そして処理部１２は、それらの状態間のエネルギー差が所定の閾値より小さくなる範囲で、磁場の強度を変更する。これにより、イジングモデルの状態が励起状態に遷移することを抑止することができる。
【0033】
〔第２の実施の形態〕
第２の実施の形態は、２値Ｎ変数で定義された最小値求解問題について系の基底状態（最適値）および基底状態の波動関数（最適値における状態変数の組み合わせ）を求めるものである。
第２の実施の形態は、２値Ｎ変数で定義された最小値求解問題について系の基底状態（最適値）および基底状態の波動関数（最適値における状態変数の組み合わせ）を求めるものである。
【0034】
図２は、第２の実施の形態のシステム構成の一例を示す図である。例えば古典コンピュータ１００に量子アニーリング型量子コンピュータ２００が接続されている。古典コンピュータ１００は、ノイマン型コンピュータとも呼ばれる。古典コンピュータ１００は、量子アニーリング型量子コンピュータ２００を制御して、イジングモデルを用いた量子アニーリングによる最適値の探索を実行させる。量子アニーリング型量子コンピュータ２００は、量子アニーリングによって組合せ最適化問題の解を探索する装置である。
【0035】
＜古典コンピュータのハードウェア構成＞
図３は、古典コンピュータのハードウェアの一例を示す図である。古典コンピュータ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。
図３は、古典コンピュータのハードウェアの一例を示す図である。古典コンピュータ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。
【0036】
メモリ１０２は、古典コンピュータ１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に利用する各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。
【0037】
バス１０９に接続されている周辺機器としては、ストレージ装置１０３、ＧＰＵ１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。
【0038】
ストレージ装置１０３は、内蔵した記録媒体に対して、電気的または磁気的にデータの書き込みおよび読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。
【0039】
ＧＰＵ１０４には、モニタ２１が接続されている。ＧＰＵ１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、有機ＥＬ（Electro Luminescence）を用いた表示装置や液晶表示装置などがある。
【0040】
入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。
【0041】
光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取り、または光ディスク２４へのデータの書き込みを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ－ＲＡＭ、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ－Ｒ（Recordable）／ＲＷ（ReWritable）などがある。
【0042】
機器接続インタフェース１０７は、古典コンピュータ１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。
【0043】
ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。
【0044】
古典コンピュータ１００は、以上のようなハードウェアによって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示した情報処理装置１０も、図３に示した古典コンピュータ１００と同様のハードウェアにより実現することができる。
【0045】
古典コンピュータ１００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。古典コンピュータ１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、古典コンピュータ１００に実行させるプログラムをストレージ装置１０３に格納しておくことができる。プロセッサ１０１は、ストレージ装置１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。また古典コンピュータ１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばプロセッサ１０１からの制御により、ストレージ装置１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。
【0046】
次に、古典コンピュータ１００が有する機能について説明する。
図４は、古典コンピュータが有する機能を示すブロック図である。古典コンピュータ１００は、記憶部１１０と処理部１２０とを有する。記憶部１１０は、例えばメモリ１０２またはストレージ装置１０３の記憶領域の一部によって実現される。処理部１２０は、例えばプロセッサ１０１に所定のプログラムを実行させることにより実現される。
図４は、古典コンピュータが有する機能を示すブロック図である。古典コンピュータ１００は、記憶部１１０と処理部１２０とを有する。記憶部１１０は、例えばメモリ１０２またはストレージ装置１０３の記憶領域の一部によって実現される。処理部１２０は、例えばプロセッサ１０１に所定のプログラムを実行させることにより実現される。
【0047】
記憶部１１０は、エネルギー情報、スピン情報、横磁場情報、非線形磁場情報、問題設定情報、ハミルトニアン情報、量子相転移情報を記憶する。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーの値のうち、最小値から小さい順に所定数のエネルギーの値を含む。また、エネルギー情報は、所定数のエネルギーの値に対応した各状態変数の値の組合せを含んでいてもよい。スピン情報は、各状態変数の値や、スピン初期化手法の指定情報を含む。横磁場情報は、系に作用する横磁場の強度の情報を含む。非線形磁場情報は、系に作用する非線形の横磁場の強度の情報を含む。問題設定情報は、例えば、求解対象の問題を示す情報、使用する最適化方法（ＳＡ法、ＱＴＳＡ（Quantum Thermal Simulated Annealing）法など）、最適化方法の実行に用いる各種パラメータ、計算終了条件などを含む。ハミルトニアン情報は、例えばエネルギー関数の重み係数（Ｗij）、バイアス係数（ｂi）、定数（Ｃ）などを含む。量子相転移情報は、量子相転移が発生する横磁場の強度などの量子相転移に関する情報を含む。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーの値のうち、最小値から小さい順に所定数のエネルギーの値を含む。また、エネルギー情報は、所定数のエネルギーの値に対応した各状態変数の値の組合せを含んでいてもよい。スピン情報は、各状態変数の値や、スピン初期化手法の指定情報を含む。横磁場情報は、系に作用する横磁場の強度の情報を含む。非線形磁場情報は、系に作用する非線形の横磁場の強度の情報を含む。問題設定情報は、例えば、求解対象の問題を示す情報、使用する最適化方法（ＳＡ法、ＱＴＳＡ（Quantum Thermal Simulated Annealing）法など）、最適化方法の実行に用いる各種パラメータ、計算終了条件などを含む。ハミルトニアン情報は、例えばエネルギー関数の重み係数（Ｗij）、バイアス係数（ｂi）、定数（Ｃ）などを含む。量子相転移情報は、量子相転移が発生する横磁場の強度などの量子相転移に関する情報を含む。
【0048】
処理部１２０は、記憶部１１０内のデータに基づいて、量子アニーリング型量子コンピュータ２００を制御して最小値求解処理を行う。例えば処理部１２０は、記憶部１１０から上記の各種情報を、処理部１２０が理解可能な形式で読み込む。また処理部１２０は、問題設定情報に基づいて、求解対象の問題を求解するためのイジングモデルのハミルトニアンを示すハミルトニアン情報を生成する。生成されるハミルトニアンには、非線形の横磁場高が含まれる。処理部１２０は、生成したハミルトニアンを記憶部１１０に格納する。
【0049】
処理部１２０は、実数時間の発展に伴って更新される横磁場（後述のＧ2）を計算する。例えば処理部１２０は、予め定義されたアニーリングスケジュールに従って、外部磁場を減少させる。その際、処理部１２０は、量子相転移の発生前後で、異なるアニーリングスケジュールを適用する。処理部１２０は、例えば量子相転移が発生する横磁場の強度の設定入力を受け、指定された横磁場の強度を境にアニーリングスケジュールを変更する。また処理部１２０は、アニーリングの過程で量子相転移の発生を検知し、量子相転移の発生を検知した後は、それ以前と異なるアニーリングスケジュールを適用することもできる。
【0050】
また処理部１２０は、虚数時間の発展を制御する。例えば処理部１２０は、横磁場を固定した状態でスピンの状態が計算されるごとに、虚数時間の時間刻み幅Δτずつ虚数時間を進める。さらに処理部１２０は、初期化されたスピン状態にノイズを付加する。例えば処理部１２０は、複数のスピンのうちの所定数のスピンをランダムに選択する。処理部１２０は、選択したスピンの確率振幅に対してノイズを付加する。処理部１２０は、スピンごとのノイズの量を、例えば乱数を用いて個別に決定する。
【0051】
なお、図４に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。また、図４に示した各要素の機能は、例えば、その要素に対応するプログラムモジュールをコンピュータに実行させることで実現することができる。
【0052】
このような古典コンピュータ１００を用いて、２値Ｎ変数で定義された最小値求解問題について、量子力学に原理を置く計算を行い、系の基底状態および励起状態の波動関数を求めることができる。ただし、前述のように、ＱＡ法を用いた解探索の時間が長期化する原因の１つとして、１次の量子相転移と呼ばれる現象の発生がある。すなわち１次の量子相転移の発生を抑制しようとすると、解探索の時間が長期化してしまう。そこで、１次の量子相転移の発生を抑止しようとすると、最小値求解問題の解探索時間が長期化する理由について、詳細に説明する。
【0053】
なお、イジングモデルの基底状態とは、イジングモデルを表すハミルトニアンのエネルギーが最も低い状態である。以下、イジングモデルの基底状態を、ハミルトニアンの基底状態と呼ぶこともある。また、ＱＡ法においてイジングモデルに印加される磁場は横磁場と呼ばれる。
【0054】
＜基底状態の算出手法の概要＞
量子アニーリングにより、古典イジングモデルを解いた場合、既存のマルコフ連鎖モンテカルロ法による求解よりも指数関数的に高速に解けることが期待されている。ただし、これは多項式時間で解けることを直ちに意味しない。一般的に２Nの計算量が必要な問題を、指数関数的に高速化した結果Ｎ1＜ＮなるＮ1を用いて２N1（「Ｎ１」の「１」は下付きの添字）の計算量になったとする。この場合において、計算速度は２N-N1（「Ｎ１」の「１」は下付きの添字）だけ加速し、確かに指数関数的な加速になるが、加速したとしても、これは多項式時間ではない。しかしながら、例え不幸にも多項式時間で解けなかったとしても、指数関数的な加速には意味がある。例えば、ある種の問題が１３８億年というふうに、宇宙開闢以来の時間を掛けても解けなかったとする。この場合に、我々はこの問題を解けないと表現する。しかし、これが不完全ながらも指数関数加速を行うことができ、１３８億倍だけ加速出来れば、計算時間は１年になる。この場合、一般には指数関数的な計算量を持つ問題を解いたことにはならないが、特定の問題を解くには十分な速度向上が得られる場合があり、実用性において技術を開発する具体的なメリットがある。
量子アニーリングにより、古典イジングモデルを解いた場合、既存のマルコフ連鎖モンテカルロ法による求解よりも指数関数的に高速に解けることが期待されている。ただし、これは多項式時間で解けることを直ちに意味しない。一般的に２Nの計算量が必要な問題を、指数関数的に高速化した結果Ｎ1＜ＮなるＮ1を用いて２N1（「Ｎ１」の「１」は下付きの添字）の計算量になったとする。この場合において、計算速度は２N-N1（「Ｎ１」の「１」は下付きの添字）だけ加速し、確かに指数関数的な加速になるが、加速したとしても、これは多項式時間ではない。しかしながら、例え不幸にも多項式時間で解けなかったとしても、指数関数的な加速には意味がある。例えば、ある種の問題が１３８億年というふうに、宇宙開闢以来の時間を掛けても解けなかったとする。この場合に、我々はこの問題を解けないと表現する。しかし、これが不完全ながらも指数関数加速を行うことができ、１３８億倍だけ加速出来れば、計算時間は１年になる。この場合、一般には指数関数的な計算量を持つ問題を解いたことにはならないが、特定の問題を解くには十分な速度向上が得られる場合があり、実用性において技術を開発する具体的なメリットがある。
【0055】
さらに、量子コンピュータの社会応用の検討も始まっており、量子コンピュータそのものを作るだけではなく、量子コンピュータをどのような社会的な課題の解決に適用していくのかという応用適用の議論も盛んになりつつある。例えば、ヘルスケア領域における病院の予約の最適化などは一つの具体的なターゲットになりえる。なぜならば、病院の予約等のスケジュールを高速で最適化できれば、病院は限られた医師等の人的リソースを最適化でき、経営効率も上がる。同時に患者側も待ち時間が減るという具体的なメリットが双方に存在する。両者が利益を有する関係になるため、量子コンピュータの応用適用先としてスケジュール最適化などの技術を開発するための具体的な動機が存在する。
【0056】
イジングモデルの基底状態をＱＡ法で求解する場合、横磁場を加えた項が初期状態として用意される。熱的散逸処理を行わない場合、ここから、求めたい２値Ｎ変数問題として定義された最適化問題のハミルトニアンの横磁場の強さを十分ゆっくり弱めていき、最終的に横磁場が０になるようにする。このとき、系の波動関数を初期状態から求解対象のイジングモデルの基底状態の波動関数に十分ゆっくり遷移させていき、励起状態に励起させないように磁場を弱めることとなる。これを断熱遷移と言う。断熱遷移の結果として、最終的に系は求解対象のハミルトニアンの基底状態に遷移する。このようにして求解対象の量子イジングハミルトニアンの基底状態を求めることができる。
【0057】
一方でＱＡにも欠点が存在する。系のハミルトニアンによっては基底状態と励起状態のエネルギーギャップが小さくなり、断熱遷移にかける時間が伸びてしまう問題がある。具体的には断熱遷移をしていく過程で、中間のハミルトニアンにおける基底状態と第一励起状態のエネルギー差が小さくなってしまう問題である。
【0058】
エネルギー差がスピンの数の関数として指数関数的に小さくなる場合を１次の量子相転移と呼ぶ。量子相転移が１次の場合、ＱＡによる断熱遷移条件を守りながら時間発展するための時間刻み幅は基底状態と励起状態のエネルギーギャップに依存する。そのため、時間刻み幅を指数関数的に小さくすることが求められる。これは量子アニーリング法により終状態に到達するための時間発展回数が指数関数的に延びることを意味する。
【0059】
他方において、エネルギー差がスピンの数の関数としてべき乗で小さくなる場合も存在する。これを２次の量子相転移と呼ぶ。２次の量子相転移が起こる場合、ＱＡの断熱遷移の速度はべき乗で小さくなる。この場合は終状態にたどり着くために必要な時間は多項式時間で増加する。
【0060】
しかし、応用上興味のある問題の多くが１次の量子相転移を伴うため、多くの実用上の問題において計算時間が指数関数的に増大する。そこで１次の量子相転移を２次の量子相転移に置き換えることが可能であるかを探る研究もなされている。Non-Stoquastic Hamiltonianと呼ばれる特殊なハミルトニアンの系において、横磁場項に工夫をすることで１次相転移を２次相転移にすることができたという報告もある（非特許文献１）。従って、通常のＱＡ法の枠組みにおいては１次の量子相転移をいかにして回避するかが計算を実行する上での中心的な課題になる。
【0061】
さらにＱＡにおいては断熱遷移条件を守ることが求められている。しかし、ＱＡ法に基づく量子コンピュータを用いた計算では、この断熱条件が破られているにも拘らず、有限の確率で基底状態を計算できている可能性があることが報告されている（非特許文献２）。量子アニーリングの実機を用いた計算においては、最適なＴＴＳ（Time To Solution）が存在するという報告であるが、物理学的な見地において、なぜ最適なＴＴＳが存在するかについては一考の余地がある。
【0062】
これらの状況を考えると、量子アニーリングにおいては、何らかの原因でエネルギー散逸（熱的散逸）が起こり、エネルギー散逸によって量子アニーリングの機構が機能している可能性があると考えられる。時間依存シュレーディンガー方程式を虚数時間で解くことにより、現象論的なエネルギー散逸効果を取り入れた計算法も提案されている（特許文献１）。
【0063】
しかしながら、エネルギー散逸効果を取り入れた計算法を用いても、基底状態を計算できない場合があり得る。それは、一度スピンが反転し、対応する古典的な状態に量子状態が陥ってしまうと、そのあとは、幾ら横磁場を弱めてもスピンは反転しないためである。つまり、巡回セールスマン問題などで、最適解以外の状態（励起状態）が虚数時間発展で得られてしまうと、励起状態から基底状態へ脱励起させるメカニズムが存在しない。
【0064】
この問題を簡単に説明する。まず、量子アニーリングハミルトニアンは通常次の（１．１）式で定義される。なお、スピン数をＮ（Ｎは自然数）とする。
【0065】
【数1】
【0066】
Ｇ1，Ｇ2は定数である。量子アニーリングではＧ1＝１－Ｇ2とする場合が多い。Ｇ1＝１－Ｇ2とする場合は、０≦Ｇ1≦１、０≦Ｇ2≦１とするが、一般には任意でよい。Ｇ1＝１－Ｇ2なる条件を課さない場合においても、境界条件としてＧ2＝０の時に、Ｇ1＝１が課される。ＨTHは最小値求解対象のハミルトニアンである。Ｔ（＾付き）は横磁場演算子である。古典イジングモデルでは２値型の変数ｚi＝｛０，１｝に対して、Ｈ（｛ｚ｝）を（１．２）式で定義する。
【0067】
【数2】
【0068】
この（１．２）式が求解対象ハミルトニアンである。ここで、Ｗij、ｂiおよび、Ｃは実数であり、それぞれ、重み係数、バイアス係数、定数項と呼ばれる。（１．２）式はそのままでは量子力学で取り扱うことができない。そのため、（１．２）式において、２値変数ｚiについて、ｚi→σi
z（σは＾付き）の置き換えをして量子化された（１．３）式のハミルトニアンＨTH（Ｈは＾付き）が導入される。
【0069】
【数3】
【0070】
ただし、このままではスピンの運動を記載する項がない。通常のシュレーディンガー方程式における運動エネルギー項がないため、横磁場項と呼ばれるスピンの運動を記述する項Ｔ（＾付き）を（１．４）式の形で導入する。
【0071】
【数4】
【0072】
つまり、（１．２）式の最小値求解をするため、ハミルトニアンを量子力学的なハミルトニアンの（１．１）式に拡張し、（１．１）式の最低エネルギー状態を求めることで、古典イジングモデルの（１．２）式の最小値求解を達成する戦略である。このような戦略を採るのは、量子論に拡張したほうが高速に最低エネルギー状態を求められると考えることができるためである。
【0073】
（１．１）式のハミルトニアンに対する、非相対論的なシュレーディンガー方程式は（１．５）式となる。
【0074】
【数5】
【0075】
量子アニーリングの理論の枠組みでは、初期状態として実験的に作成可能な状態を初期状態にすることになる。具体的には（１．６）式となる。
【0076】
【数6】
【0077】
ここで、｜φi（ｔ＝０）＞はｉ＝１，２，・・・，Ｎに対して定義され、１スピンの波動関数である。量子アニーリングでは、（１．７）式が初期状態とされる。
【0078】
【数7】
【0079】
これは、（１．７）式がＪi＞０において、横磁場演算子Ｔ（＾付き）の基底状態になっているからである。（１．７）式の初期状態から系の波動関数を（１．５）式に従い時間発展する処理が行われる。つまり、Ｇ2（ｔ＝０）＝１、Ｇ1（ｔ＝０）＝０と設定される。励起状態に励起させないように、Ｇ2を十分ゆっくり小さくし、Ｇ2（ｔ≫１）＝０となるようにする。このとき、Ｇ1（ｔ≫１）＝１となるので、励起状態に遷移しないのであれば、最終的には求解ハミルトニアンＨTH（Ｈは＾付き）の基底状態が得られる。このとき、励起状態に遷移しない条件を断熱遷移条件と呼ぶ。
【0080】
なお、演算子は以下のように定義される。
【0081】
【数8】
【0082】
【数9】
【0083】
【数10】
【0084】
【数11】
【0085】
（１．８）式は通常用いられるスピン演算子の定義と異なっているが、これはスピン演算子の固有値を１および－１ではなく、１および０にするためである。｜αi＞と｜βi＞には直交関係として（１．１２）式が成立する。ここでδijはクロネッカーのデルタである。
【0086】
【数12】
【0087】
次に縦磁場項σi
zについての固有方程式を以下に示す。
【0088】
【数13】
【0089】
【数14】
【0090】
横磁場項σi
x（σは＾付き）については以下のようになる。
【0091】
【数15】
【0092】
【数16】
【0093】
留意すべきはσi
z（σは＾付き）演算子の定義が（１．９）式であり、通常用いられるパウリのスピン行列と若干異なっていることである。
前述の特許文献１で導入したＱＴＳＡ法は（１．５）式のシュレーディンガー方程式を数値的に解く場合に、時間を実数から虚数にして解く方法である。時刻ｔにおける波動関数は（１．１７）で一般に与えられる。
前述の特許文献１で導入したＱＴＳＡ法は（１．５）式のシュレーディンガー方程式を数値的に解く場合に、時間を実数から虚数にして解く方法である。時刻ｔにおける波動関数は（１．１７）で一般に与えられる。
【0094】
【数17】
【0095】
ここで、｜ψn＞はＨTH（Ｈは＾付き）の固有状態とする。すなわち、固有状態が（１．１８）である。
【0096】
【数18】
【0097】
ＱＴＳＡ法では、ここで時間ｔが（１．１９）式を用いて虚数に変換される。
【0098】
【数19】
【0099】
τは実数であるが、虚数時間と呼ばれる。（１．１９）式を用いて（１．１７）式を変形すると次の（１．２０）式が得られる。
【0100】
【数20】
【0101】
（１．２０）式において、虚数時間τ→∞の極限を取ると、基底状態のみの振幅が残り、（１．２１）式のようになる。
【0102】
【数21】
【0103】
ここで、（１．２１）式も最終的には０になってしまうが、これは時間を虚数にすることでユニタリ保存が破れるためである。しかし、（１．２０）式において、強制的にユニタリ保存するように、波動関数の規格化条件を守るように再規格化しながら計算すると、励起状態の情報のみが消えて、最終的に基底状態の情報のみが残る。これが虚数時間発展法の概要である。具体的な手続きは特許文献１に記載された通りである。
【0104】
しかし、虚数時間発展の前提になる完全系展開では（１．１７）式に「基底状態が含まれている」ことを前提とする。虚数時間発展の途中で基底状態の情報が失われてしまえば、励起状態に収束してしまう。そのため、一度励起状態に収束してしまうと、虚数時間発展を続けても基底状態に落ちることはない。これは、虚数時間発展の前提となる基底状態の情報が消えてしまうのが原因である。
【0105】
図５は、虚数時間発展法を用いた熱的散逸の模式図である。図５には、横軸にλ、縦軸にハミルトニアンによって求められる系のエネルギーをとったグラフ３０が示されている。図５の例では、実数時間発展に伴い横磁場の強さに対応する定数λが０から１に遷移するときの、λに応じた基底状態におけるエネルギーを実線３１で示している。また、λに応じた第１励起エネルギー状態を点線３２で示している。λが０のときに横磁場の強さが最大であり、λが１のときに横磁場の強さが最小である。
【0106】
ここで図５中のλは１－Ｇ2であるとする。ある与えられたλに対して、虚数時間発展を行い、与えられたλにおける基底状態を求めるものとする。そして、λの値を少し増やしてλ＋Δλとすることを想定する。このとき、λ＋Δλにおける波動関数の初期値として、前段のλの値における虚数時間発展の結果得られた波動関数が用いられる。この手続きの繰り返しにより、最終的にλ＝１における波動関数を求めることができる。
【0107】
しかし、１次の量子相転移が起こる状況では図５における点Ｐ付近のように、複数の量子状態がエネルギー的に近接している。そのため、Δλを十分小さく取らないと、励起状態に遷移してしまうのである。つまり、本来は点Ｑから点Ｒに虚数時間発展の結果遷移することを期待するが、一旦点Ｑに遷移してしまうと、これは励起状態である。そのため、遷移した時点で基底状態の情報が失われてしまう。結果として、虚数時間発展を続けても点Ｒに遷移しない。そして、λ＝１まで虚数時間発展を続けても、もはや基底状態に到達することはない。
【0108】
実際には量子ビット数を増やすと、エネルギー的に近接した量子状態が非常に多くなる。そして、エネルギー差も非常に小さくなってしまう。そのため、系にほんの少しノイズが乗るだけで、最終的に得られるλ＝１の状態は変わってしまう。このとき、スピン配列はかなり変わった状態が得られる。エネルギー差は問題によっては（１．２２）式のように指数関数的に小さくなることが知られているが、中間状態のほんのわずかなエネルギー差により、結果として得られる終状態が大きく変わってしまう。
【0109】
【数22】
【0110】
計算や実験は必ず有限精度で行われることになるため、有限精度の測定では基底状態と励起状態が区別できなくなる。例えば、シミュレーションでは倍精度実数を使うことが多いが、倍精度実数は１６桁しか扱うことができない。もし、基底状態と励起状態がほぼ縮退しており、そのエネルギー差が１７桁目以降に出るような状態であれば、１６桁の精度の計算ではもはや区別することができない。数値的には縮退した状態になってしまう。
【0111】
結果として、数値的に縮退した状態と見なされてしまうと、どの状態に遷移するかは完全にランダムに決まってしまう。波動関数の時間発展時に現れたわずかな誤差で遷移先が決定されてしまう。そして、問題はその遷移先をもはや制御できないことである。
【0112】
これは、自発的対称性の破れが起こっていると考えればよい。複数の実質的に縮退したエネルギー準位構造があり、量子アニーリング時には、どのエネルギー準位にも等確率に遷移できてしまう。行き先を決めるのは計算時に蓄積した誤差や外界からのノイズなどほんの少しの初期状態の違いによって、最終的に行き着く巨視的な量子状態が変わってしまう。偶発的な要因で系の運命は決まってしまう。
【0113】
では３２桁や６４桁の多倍精度を使えば解決するのだろうか。（１．２２）式が正しいのであれば、つまり、基底状態と励起状態のエネルギー差がスピン数Ｎの関数として指数関数的に小さくなるのであれば、スピン数を増やしたときにはエネルギー差を区別するために必要な桁数はスピン数Ｎに比例して増えていくことになる。そのため、Ｎ≫１の場合、エネルギー差を記述するために用いられる桁数ＮdはＮd∝Ｎであり、たとえ、３２桁や６４桁の多倍長演算を用いても、倍精度計算時と同じ問題が起こってしまう。
【0114】
加えて、もう一つ大きな問題がある。計算時間の問題である。量子アニーリング時の計算時間（虚数時間発展の計算回数）ＮIは、ある与えられたＧ2の値での虚数時間発展に要する計算回数ｎI（Ｇ2）の関数として、（１．２３）式で与えられる。
【0115】
【数23】
【0116】
ここで、Ｇ20は１次の量子相転移が起こるＧ2の値である。虚数時間発展における計算時間の見積もりは、（１．２０）式より、第一励起状態が十分小さくなる時間として与えられるため、（１．２４）式で与えられる。
【0117】
【数24】
【0118】
このため、計算時間τ1は（１．２５）式で見積もることができる。
【0119】
【数25】
【0120】
ΔＥ1は（１．２６）式の通りである。
【0121】
【数26】
【0122】
ここで１次の量子相転移が起こるのであれば、（１．２２）式でΔＥ1が与えられるため、（１．２７）式に示す程度の計算時間がかかる。
【0123】
【数27】
【0124】
従って、全体の計算時間ＮIは（１．２８）式で表される。
【0125】
【数28】
【0126】
これより、１次の量子相転移が起こると、計算時間はスピン数Ｎの関数として指数関数的に長くなってしまう。しかも、（１．２５）式はエネルギーと時間の不確定性原理と呼ばれるものである。このため、（１．２５）式は計算技術上の問題ではなく、物理学としての原理的な問題である。
【0127】
この問題の根源は、エネルギー差が（１．２２）式のように指数関数的に小さくなるという所から来ており、実質上エネルギー準位が縮退してしまうことが問題の本質である。この場合、エネルギー準位が実質的に縮退してしまうため、ユーザが計算の最終状態を制御するのが困難になる。この問題は、１回の量子アニーリング計算で基底状態を得ようとする場合において発生する。
【0128】
このような状況にあっても、確率的に基底状態を求めることができればよいのであれば、短時間で基底状態を求めることは可能である。
そこで導入するのが量子試行の概念である。ここで、虚数時間発展で実際に４都市の巡回セールスマン問題（ＴＳＰ：Traveling Salesman Problem）を解く場合を例にとり、量子試行の概念を簡単に説明する。
そこで導入するのが量子試行の概念である。ここで、虚数時間発展で実際に４都市の巡回セールスマン問題（ＴＳＰ：Traveling Salesman Problem）を解く場合を例にとり、量子試行の概念を簡単に説明する。
【0129】
図６は、４都市のＴＳＰにおける都市の配置例を示す図である。図６の例では、１番目の都市の位置が（０，０）、２番目の都市の位置が（３，０）、３番目の都市の位置が（３，１）、４番目の都市の位置が（１，１）である。この問題では経路が３つしかなく、最小の経路から順に７．４１４、７．８１２、１０．３９８であることが分かっている。
【0130】
初期状態として、（１．２９）式のように波動関数の振幅にほんの少しだけノイズδiを入れる。ただし、ノイズを付加すると波動関数のノルムが１から変化するため、ノイズ付加後に、波動関数のノルムが１になるように再規格化する。
【0131】
【数29】
【0132】
（１．２９）式のノイズはランダムな値が設定される。ただし、虚数時間発展時の収束判定は計算に用いる桁数よりも小さく取っておく。例えば、倍精度浮動小数点を使用する場合、エネルギーの収束判定はΔＥ≦１０-10などとしておく。このようにすると、初期に与えたノイズの影響が残ったままになる。このようにして、系にわずかにランダムにノイズを与えておき、量子アニーリング計算を行う。計算結果として一つの量子状態が最終的に得られる。この一連の操作を量子試行と呼ぶことにする。
【0133】
量子試行を２００回ほど実行して得られた結果を図７に示す。ここで、Ｇ2を離散的な点として与え、等間隔に分点を定義する。そして、計算に用いる分点の数は１０２４である。
【0134】
図７は、初期状態にランダムに付加したノイズに応じた計算結果の違いを示す図である。図７に示すグラフ３３は横軸がＧ2であり、縦軸が全エネルギーである。グラフ３３に示す３本の実線それぞれが、異なる分点数の量子試行によるＧ2の値ごとのエネルギー値を示している。グラフ３３に示すように、分点の数によっては結果が７．４１４の基底状態が求まる場合もあれば、７．８１２や１０．３９８の励起状態が求まる場合もある。
【0135】
図７に示すように、初期状態として与えたわずかなノイズにより最終的に得られるエネルギー状態も異なってしまうが、有限の確率で基底状態に到達する。これは、問題を一種の散乱過程として捉えなおして解いていると考えることができる。
【0136】
しかしながら、量子試行を用いた場合の課題も存在する。量子試行を用いた計算によって、量子アニーリングのためのアニーリングスケジュールは大幅に緩和可能であるし、量子試行によって有限の確率で基底状態も計算ができる。しかし、スピン数Ｎが増えると、基底状態が得られる確率が下がってしまい、基底状態に到達することが難しくなる。
【0137】
このような背景から、量子試行を実行するときに、基底状態に到達する確率を高めるための技術開発が望まれている。
そこで古典コンピュータ１００では、量子試行を用いた場合における基底状態に到達する確率を高めるための技術が適用される。
そこで古典コンピュータ１００では、量子試行を用いた場合における基底状態に到達する確率を高めるための技術が適用される。
【0138】
＜基底状態に到達する確率を高めるための技術の概要＞
古典コンピュータ１００は、２値型Ｎ変数で定義されるイジングモデルの最適解を量子力学に基づく変分原理で求める。そして古典コンピュータ１００は、提案の技術を適用することで、ＱＴＳＡ法で計算をした場合、必ずしも基底状態が求まらず、場合によっては基底状態にたどり着かず、励起状態に到達してしまう問題を解決し、高確率で基底状態にたどり着くようにする。
古典コンピュータ１００は、２値型Ｎ変数で定義されるイジングモデルの最適解を量子力学に基づく変分原理で求める。そして古典コンピュータ１００は、提案の技術を適用することで、ＱＴＳＡ法で計算をした場合、必ずしも基底状態が求まらず、場合によっては基底状態にたどり着かず、励起状態に到達してしまう問題を解決し、高確率で基底状態にたどり着くようにする。
【0139】
古典コンピュータ１００は、基底状態に到達する確率を高めるために行う処理は、例えば以下のような処理である。
１．古典コンピュータ１００は、量子アニーリングハミルトニアンを量子力学の枠組みで解く。すなわち古典コンピュータ１００は、熱散逸の効果を取り入れ、量子力学的な時間発展演算子における時間を複素数に拡張する。そして古典コンピュータ１００は、虚数時間発展法に基づき、熱散逸を実現する。熱散逸機構の導入により、断熱遷移条件を厳密に守らずとも基底状態へ緩和できる。
１．古典コンピュータ１００は、量子アニーリングハミルトニアンを量子力学の枠組みで解く。すなわち古典コンピュータ１００は、熱散逸の効果を取り入れ、量子力学的な時間発展演算子における時間を複素数に拡張する。そして古典コンピュータ１００は、虚数時間発展法に基づき、熱散逸を実現する。熱散逸機構の導入により、断熱遷移条件を厳密に守らずとも基底状態へ緩和できる。
【0140】
２．古典コンピュータ１００は、ＱＴＳＡ法を虚数時間における散乱過程として捉えなおし、初期状態に多様性を与える。このことで、初期状態と終状態のペアができる。これを量子試行と定義し、古典コンピュータ１００は、初期状態を変えた多数回の量子試行によって対象ハミルトニアンの基底状態を含む励起状態を求める。
【0141】
３．古典コンピュータ１００は、波動関数の初期値に多様性を持たせるため、通常は邪魔者扱いされるノイズを積極的に利用し、ノイズにより系の波動関数に多様性を持たせる。その結果として、系の対称性が破壊され、量子アニーリングの計算における計算量の致命的な増大が緩和される。
【0142】
４．古典コンピュータ１００は、スピン数Ｎが大きくなると、基底状態へ到達する確率が低くなるため、非線形の外部磁場を入れる。これにより、多スピン励起を誘発し、実質上のエネルギー縮退が起こっている量子状態の数を減らすことができる。そのため、終状態の候補数が少なくなる。結果として基底状態への遷移確率を高めることができる。
【0143】
以下、古典コンピュータ１００が実行する処理について詳細に説明する。
＜イジングモデル問題のハミルトニアンについて＞
まず、最小値求解問題であるイジング型問題を説明する。イジング型問題は元々物理学領域において磁性体の研究に使われたモデルであり、Ｎ個のスピンがあるときの全エネルギーは（１．２）式で記述される。
＜イジングモデル問題のハミルトニアンについて＞
まず、最小値求解問題であるイジング型問題を説明する。イジング型問題は元々物理学領域において磁性体の研究に使われたモデルであり、Ｎ個のスピンがあるときの全エネルギーは（１．２）式で記述される。
【0144】
イジング型問題は、（１．２）式で与えられるハミルトニアンの最小値を求めるものである。しかし、最小値を求めるための計算量は一般に２N回程度になることが知られている。つまり、最適解を求めるための計算時間は多項式時間で抑えられず、指数関数的に増大する。そのため、（１．２）式を総当たりで解くようなことは通常せず、メトロポリス法に基礎を置く計算を実行して最小値を求めることが多い。ＳＡ法やレプリカ交換法などがよく使用されるが、メトロポリス法に基礎を置く計算法では、最小値にたどり着くまでの計算量が厳密に求められているわけではない。そしてメトロポリス法に基礎を置く計算法は、経験的にはスピンの数Ｎに対して計算量が指数関数的に増大することが知られている。
【0145】
古典的なハミルトニアンを確率過程に基礎を置く計算よりも高速化するため、量子論を用いて計算を行う枠組みを考える。ただし、ハミルトニアンの（１．１）式から少し変えて、（２．１）式で定義される量子力学的なハミルトニアンを用いる。
【0146】
【数30】
【0147】
（２．１）式の（１．１）式との違いは横磁場に加えて、付加的な外部磁場項が加えられていることである。ここでは、Ｖadd（Ｖは＾付き）として非線形の外部磁場を想定する。具体的な例として（２．２）式を考える。
【0148】
【数31】
【0149】
ここで、Ｇ3については以下の境界条件（２．３）および（２．４）を満たすものとする。
【0150】
【数32】
【0151】
【数33】
【0152】
この条件を満たすものであれば、Ｇ3としてさまざまな取り方が考えられる。例えば、一つの取り方としてはｎ＞０として（２．５）式のように取ることができる。
【0153】
【数34】
【0154】
この条件は、Ｇ2→１において、方程式が自明解を持つようにするためである。Ｇ2→０における条件は、Ｇ2→０でハミルトニアンが求めるべきポテンシャルと一致するようにするためである。なお、Ｇ1については、量子アニーリングでよく使用される型を用いることができ、（２．６）式で表される。
【0155】
【数35】
【0156】
ＨTH（Ｈは＾付き）および、Ｔ（＾付き）の定義はそれぞれ（１．３）式、（１．４）式とする。演算子の定義は（１．８）～（１．９）とする。ブラとケットベクトルの定義と演算の定義は（１．１０）～（１．１６）とする。非相対論的なシュレーディンガー方程式は（１．５）とする。波動関数は１スピンの波動関数の直積形として（１．１６）式とする。
【0157】
次にｉ番目の量子ビットを表す１粒子波動関数のケット表示を（２．７）式で定義する。
【0158】
【数36】
【0159】
Ｃiα（ｔ）は、ｉ番目のスピンの｜αi〉の状態の確率振幅と位相を表す複素数である。Ｃiβ（ｔ）は、ｉ番目のスピンの｜βi〉の状態の確率振幅と位相を表す複素数である。１粒子波動関数のブラ表示を次の（２．８）式で定義する。
【0160】
【数37】
【0161】
ここで、１粒子波動関数は規格化されているとして、（２．９）式が成り立つ。
【0162】
【数38】
【0163】
次によく使用する物理量を以下のように定義する。
【0164】
【数39】
【0165】
【数40】
【0166】
Ｎ粒子系の波動関数を以下の様に定義する。
【0167】
【数41】
【0168】
なお、ｍ番目のスピンがαスピンとして観測される確率をＰm,α、βスピンとして観測される確率をＰm,βとすると、それぞれ以下の式で表される。
【0169】
【数42】
【0170】
【数43】
【0171】
以上がイジングモデル問題のハミルトニアンの説明である。
＜虚数時間発展法による基礎方程式＞
次に虚数時間発展法による基礎方程式について説明する。
＜虚数時間発展法による基礎方程式＞
次に虚数時間発展法による基礎方程式について説明する。
【0172】
今、ハミルトニアン（２．１）式で定義される時間依存シュレーディンガー方程式（１．５）式を解くことを考える。
次に虚数時間発展法（ＩＴＰ：Imaginary Time Propagation）を量子イジング系に適用した場合の支配方程式を導出する。時間発展演算子として１次のＩＴＰを用いる。１次のＩＴＰの演算子は（２．１５）式で与えられる。
次に虚数時間発展法（ＩＴＰ：Imaginary Time Propagation）を量子イジング系に適用した場合の支配方程式を導出する。時間発展演算子として１次のＩＴＰを用いる。１次のＩＴＰの演算子は（２．１５）式で与えられる。
【0173】
【数44】
【0174】
（２．１５）式をΔｔに関して１次の項まで展開して整理すると（２．１６）式となる。
【0175】
【数45】
【0176】
次に波動関数の時間発展を記述する方程式を導く。波動関数の時間発展は（２．１７）式で与えられる。
【0177】
【数46】
【0178】
今、ｍ番目の粒子の波動関数だけ除いたＮ－１スピン系の波動関数を｜Ψm（τ）＞（ｍはオーバーライン付き）と表記すると、この波動関数は次の（２．１８）式のようになる。
【0179】
【数47】
【0180】
（２．１８）式の両辺に＜αm｜×＜Ψm（τ0）（×はテンソル積（丸の中に×）、Ψの添字ｍにはオーバーライン付き）を作用させる。その結果、左辺は（２．１９）式、（２．２０）式となる。
【0181】
【数48】
【0182】
【数49】
【0183】
ここで、Ｍi
*はスピンの変化のしにくさを表す量であり、（２．２１）式で表される。
【0184】
【数50】
【0185】
ここで（２．２２）式を定義する。
【0186】
【数51】
【0187】
Ωm（τ）をΔτの項まで整理すると（２．２３）式となる。
【0188】
【数52】
【0189】
次に（２．１８）式の両辺に＜αm｜×＜Ψm（τ0）（×はテンソル積（丸の中に×））、Ψの添字ｍにはオーバーライン付き）および＜βm｜×＜Ψm（τ0）（×はテンソル積（丸の中に×））、Ψの添字ｍにはオーバーライン付き）を作用させた結果の右辺を求め、ｍ番目（ｍ＝１，２，・・・，Ｎ）のスピンの従う支配方程式を求めると（２．２４）式が求まる。
【0190】
【数53】
【0191】
ここで、行列要素は以下の（２．２５）式～（２．２７）式で与えられる。
【0192】
【数54】
【0193】
【数55】
【0194】
【数56】
【0195】
（２．２５）式～（２．２７）式中のパラメータを以下に与える。ＥTHは求解対象のハミルトニアンのエネルギー項であり、（２．２８）式～（２．３０）式で定義される。
【0196】
【数57】
【0197】
【数58】
【0198】
【数59】
【0199】
ここでｈz,mは平均場に相当する量であり、（２．３１）式および（２．３２）式で与えられる。
【0200】
【数60】
【0201】
【数61】
【0202】
横磁場項Ｔ（＾付き）起因のエネルギーは式（２．３３）の通りである。
【0203】
【数62】
【0204】
外部磁場として与える高次の横磁場項の（２．２）式のハミルトニアンのエネルギーは次の（２．３４）式で与えられる。
【0205】
【数63】
【0206】
Ｔm,ααは横磁場項の期待値であり、＜αm｜×＜Ψm（τ0）｜Ｔ｜Ψ＞（×はテンソル積（丸の中に×）、Ψmのｍはオーバーライン付き、Ｔは＾付き）からＣmα（τ）の係数を抜き出したものである。Ｔm,αβは、＜αm｜×＜Ψm（τ0）｜Ｔ｜Ψ＞（×はテンソル積（丸の中に×）、Ψmのｍはオーバーライン付き、Ｔは＾付き）からＣmβ（τ）の係数を抜き出したものである。そして、Ｔm,ββは、＜αm｜×＜Ψm（τ0）｜Ｔ｜Ψ＞（×はテンソル積（丸の中に×）、Ψmのｍはオーバーライン付き、Ｔは＾付き）からＣmβ（τ）の係数を抜き出したものである。
【0207】
Ｔm,αα
(2)、Ｔm,ββ
(2)および、Ｔm,αβ
(2)は次の（２．３５）式、（２．３６）式および、（２．３７）式で与えられる。
【0208】
【数64】
【0209】
【数65】
【0210】
【数66】
【0211】
（２．２４）式に従い、虚数時間発展を行うことになるが、このときユニタリティ（確率の保存）は破れる。そのため、時間発展計算後、波動関数を１にするように再規格化（Renormalization）が行われる。
【0212】
次にΩm（τ）について説明する。ｍ番目のスピンに対してαスピンが観測される確率をＰm,α、βスピンが観測される確率をＰm,βとする。Ｐm,αとＰm,βとは（２．３８）式、（２．３９）式で表される。
【0213】
【数67】
【0214】
【数68】
【0215】
時間発展毎に再規格化をするため、ユニタリティ（確率の保存）は強制的に保存されており、（２．４０）が満たされる。
【0216】
【数69】
【0217】
従って、観測確率の和を虚数時間τについて微分すると（２．４１）式となる。
【0218】
【数70】
【0219】
従って、Ｍm
*は実数部を持たない純虚数である。ここで、Ｍm
*はＭmの複素共役とする。
ここで、時間発展の支配方程式（２．２４）式において、波動関数の初期値としてＣmα（τ＝０）とＣmβ（τ＝０）を共に実数に選んだ場合を考える。（２．２４）式の右辺の行列要素はすべて実数である。従って、ＣmαとＣmβが実数であればＣmα（τ＋Δτ） Ωm（τ）は実数である。一方で、Ωm（τ）はＭm *から計算されるが、Ｍm *の定義式（２．２１）式から時刻τにおけるＭm *も実数になる。ところが、（２．４１）式からＭm *は純虚数になるため、この場合はＭm *＝０になる。従って、Ωm（τ）＝１になる。Ｃmα（τ＋Δτ）も実数になる。Ｃmβ（τ＋Δτ）もまったく同様の議論で実数になる。
ここで、時間発展の支配方程式（２．２４）式において、波動関数の初期値としてＣmα（τ＝０）とＣmβ（τ＝０）を共に実数に選んだ場合を考える。（２．２４）式の右辺の行列要素はすべて実数である。従って、ＣmαとＣmβが実数であればＣmα（τ＋Δτ） Ωm（τ）は実数である。一方で、Ωm（τ）はＭm *から計算されるが、Ｍm *の定義式（２．２１）式から時刻τにおけるＭm *も実数になる。ところが、（２．４１）式からＭm *は純虚数になるため、この場合はＭm *＝０になる。従って、Ωm（τ）＝１になる。Ｃmα（τ＋Δτ）も実数になる。Ｃmβ（τ＋Δτ）もまったく同様の議論で実数になる。
【0220】
次に波動関数の初期値Ｃmα（τ＝０）、Ｃmβ（τ＝０）を複素数として決めた場合を考察する。この場合、Ｍm
*は複素数として与えられる。この場合、一般にはΩm（τ）≠１であり、複素数になる。そして、Ｍm
*はｍ番目のスピンから決まる量である。従って、Ωm（τ）は着目するｍ番目以外のすべてのスピンの情報から決まる。従って、Ｃmα（τ＋Δτ）は直前のｉ≠ｍを満たす全ての他のスピンの波動関数Ｃiα（τ）とＣiβ（τ）の影響を受ける。すなわち、時刻τ＝０の波動関数の情報に位相情報を組み入れ、複素数にした場合、虚数部分を介して情報をやり取りする。絶対値の２乗で定義される観測確率には影響を与えない位相部分がスピン間の相互作用にも重要な役割を果たしている。そして、この相互作用は他の全てのスピンと結合し、位相部分を介して他のすべてのスピンともつれ合う。
【0221】
＜虚数時間発展の物理的な意味と量子力学的なシミュレーテッドアニーリング＞
虚数時間発展の枠組みにおいては、（１．１９）式のように、実数で定義された時間を虚数に変換する。このとき、量子力学的な時間発展演算子を波動関数に作用させると、τ≫１の場合には（１．２０）式のようになることは既に示した。（１．２０）式とボルツマン分布である（２．４２）式が対応すると考える。
虚数時間発展の枠組みにおいては、（１．１９）式のように、実数で定義された時間を虚数に変換する。このとき、量子力学的な時間発展演算子を波動関数に作用させると、τ≫１の場合には（１．２０）式のようになることは既に示した。（１．２０）式とボルツマン分布である（２．４２）式が対応すると考える。
【0222】
【数71】
【0223】
このときの対応関係から（２．４３）式となる。
【0224】
【数72】
【0225】
ここから、τ→∞の場合、Ｔ→０に対応することがわかる。この意味で、虚数時間発展は、量子力学的に拡張された一種のシミュレーテッドアニーリングと見なすことができる。このような解探索手法がＱＴＳＡ法である。
【0226】
＜量子アニーリングスケジュールについて＞
量子アニーリングをするときには、アニーリングスケジュールが重要と言われる。虚数時間発展法においては、Ｇ2を変える実時間発展をして、次に、エネルギー散逸を含めた虚数時間発展を行い、与えられたＧ2における最低エネルギー状態を求める過程を繰り返す。ここでＧ2を変える実時間発展のスケジュールは重要である。例えば時間発展のスケジュールはべき乗型にすることができる。すなわち、ここでは（２．４４）とする。
量子アニーリングをするときには、アニーリングスケジュールが重要と言われる。虚数時間発展法においては、Ｇ2を変える実時間発展をして、次に、エネルギー散逸を含めた虚数時間発展を行い、与えられたＧ2における最低エネルギー状態を求める過程を繰り返す。ここでＧ2を変える実時間発展のスケジュールは重要である。例えば時間発展のスケジュールはべき乗型にすることができる。すなわち、ここでは（２．４４）とする。
【0227】
【数73】
【0228】
例えばｃ（ｃは実数）はＮ-1に比例した値とする。虚数時間発展法では断熱遷移を要求しないため、ｃがＮ-1と比例関係を有していなくてもよい。
＜非線形の外部磁場の有用性＞
次にＱＡ法における非線形の外部磁場の有用性について説明する。まず、検証する課題として、Ｎスピン系のイジングモデルとして定義された巡回セールスマン問題（ＴＳＰ）を考える。
＜非線形の外部磁場の有用性＞
次にＱＡ法における非線形の外部磁場の有用性について説明する。まず、検証する課題として、Ｎスピン系のイジングモデルとして定義された巡回セールスマン問題（ＴＳＰ）を考える。
【0229】
図８は、８都市ＴＳＰの一例を示す図である。図８には８都市ＴＳＰの問題を定義する問題データ４１が示されている。問題データ４１には、８つの都市それぞれの都市番号（ｃｉｔｙ Ｎｏ．）に対応付けて、その都市の位置を示すｘ，ｙの座標が設定されている。このような８都市のＴＳＰは、６４量子ビットの系で表現される。この問題の最適解の一つの都市の訪問順序は都市１→２→３→４→５→６→７→８である。そのときの最適値は－２．５０２３である。
【0230】
図８のＴＳＰ問題を解く場合、例えば（１．１）式で定義される量子アニーリングハミルトニアンに対して、非線形の外部磁場を加えた量子アニーリングハミルトニアン（２．１）式を定義することができる。（１．１）式と（２．１）式に対して、それぞれ量子試行を４００回ずつ行い、基底状態が得られる確率が向上するかどうか検証する。
【0231】
図９は、８都市ＴＳＰの解探索結果の一例を示す図である。グラフ４２は、（１．１）式のように１次の横磁場項だけを含んだ量子アニーリングハミルトニアンに対して４００回の量子試行を行った結果を示している。グラフ４３は、（２．１）式のように２次の横磁場項を含んだ量子アニーリングハミルトニアンに対して４００回の量子試行を行った結果を示している。グラフ４２，４３は、４００回の量子試行の結果得られたエネルギーのヒストグラムを示している。グラフ４２，４３の横軸がエネルギーであり、縦軸がそのエネルギーが得られた量子試行の回数である。
【0232】
図９に示すように、非線形の横磁場項を含んだ（２．１）式を用いて量子アニーリングを行った方が、エネルギーが最小値となる割合が高い。すなわち非線形の横磁場項を含めることで基底状態が得られる確率が向上することが分かる。また、グラフ４２では、終状態として４９種の状態が得られている（得られたエネルギーが４９パターン）。グラフ４３では、終状態として２６種の状態が得られている（得られたエネルギーが２６パターン）。従って、非線形の横磁場項を含んだ（２．１）式による量子アニーリングの方が、終状態として得られる状態数が少ないことが分かる。
【0233】
＜量子相転移の発生を伴う場合のエネルギーの推移について＞
次に量子相転移の発生を伴う場合のエネルギーの推移について説明する。
図１０は、Ｇ2の関数として得られた系の全エネルギー値を示す図である。グラフ４４は、量子相転移が発生した場合のＧ2に応じたエネルギー値を示している。グラフ４５は、グラフ４４において、量子相転移が発生した後のＧ2についてエネルギーを示した図である。グラフ４５はグラフ４４の拡大図である。グラフ４４，４５の横軸がＧ2であり、縦軸がエネルギーである。
次に量子相転移の発生を伴う場合のエネルギーの推移について説明する。
図１０は、Ｇ2の関数として得られた系の全エネルギー値を示す図である。グラフ４４は、量子相転移が発生した場合のＧ2に応じたエネルギー値を示している。グラフ４５は、グラフ４４において、量子相転移が発生した後のＧ2についてエネルギーを示した図である。グラフ４５はグラフ４４の拡大図である。グラフ４４，４５の横軸がＧ2であり、縦軸がエネルギーである。
【0234】
Ｇ2は一種の断熱変数となっている。量子試行を行うと、Ｇ2＝０．７４４付近で「量子相転移」と呼ばれる現象が起こる。この値は一定ではないが、グラフ４４ではＧ2＝０．７４３±０．００１程度の非常に狭い範囲で起こっていることが分かる。そして、量子相転移が起こった後、Ｇ2の値を減らしていくと、Ｇ2＝０の極限で求解ハミルトニアンの固有状態（古典イジングモデルの状態に対応する）に収束するが、その値は基底状態と励起状態を含む状態に収束している。
【0235】
重要なことは、一旦量子相転移が起こると、もはや別の量子状態に遷移することはないことである。系の運命は量子相転移が起こるＧ2＝０．７４３±０．００１で決まってしまう。逆に、量子相転移が起こった後は、真面目にアニーリングをすることを要しない。すなわち量子相転移が起こった後はＧ2の値をかなり粗く取っても、最終結果はもはや変わらない。
【0236】
量子アニーリングのスケジュールは例えばべき乗で定義して減少させるが、このアニーリングスケジュールを厳密に守ることまでは求められない。量子相転移が起こるまでは慎重に計算すべきだが、量子相転移が起きた後は系の運命は既に決まっている。従って、量子相転移が起こった後の計算は、量子相転移が起こるＧ2の値をＧ20としたとき、例えば（３．１）式のように線形補完で決めることができる。
【0237】
【数74】
【0238】
ここで、ｍは分割数を示す自然数である。ｎは量子相転移後の何回目のＧ2の更新かを示す自然数であり、ｎ＝１，２，・・・，ｍとする。なお、極論すれば、１次量子相転移が起こった後は系の運命は決しているため、ｍ＝２にしてもよい。すなわち量子アニーリングのスケジュールを真面目に守れば、Ｇ2が小さいときほど大きい計算量が要求されることになるが、ＱＴＳＡ法の枠組みではそのような要求はない。そのため、量子相転移が発生した後は（３．１）式のようにして計算量を減らすことができる。
【0239】
この結果を量子アニーリング型量子コンピュータの実機について当てはめると、アニーリングスケジュールについては以下のパラメータを決めればよいことが分かる。横磁場である外部磁場の掃引開始時刻ｔstart、掃引終了時刻ｔend、そして、磁場のアニーリングスケジュールを変更する時刻ｔqpt、および、アニーリングスケジュール関数ｓ（ｔ）である。ただし、ｔstartにおいてＧ2＝１．０、ｔendにおいてＧ2＝０、ｔqptにおいてＧ2＝Ｇ20とする。Ｇ20は解くべき問題に依存して変わり、量子相転移が起こり、スピンが一斉に反転する横磁場項の強さＧ2がＧ20である。グラフ４４の例では、エネルギーが急激に下がる部分のＧ2の値がＧ20である。区間［ｔstart，ｔqpt］では、Ｇ2の値はＧ2（ｔ）＝ｓ（ｔ）に従い変化し、１．０＝ｓ（ｔstart）、Ｇ20＝ｓ（ｔend）を満たす。区間（ｔqpt，ｔend］では、別の高速なアニーリングスケジュールＧ2（ｔ）＝ｓ1（ｔ）を用いてよい。
【0240】
図１１は、外部磁場の時間変化の一例を示す図である。グラフ４６は、横軸にシミュレーション上の時刻と外部磁場Ｇ2との関係を示している。グラフ４６の横軸が時刻であり、縦軸が外部磁場Ｇ2である。時刻は、例えば外部磁場Ｇ2を更新する単位時間だけ進行するものとする。グラフ４６では、量子相転移が発生する時刻（図１０にも示したグラフ４４におけるエネルギーが急激に下がる時刻）を境に、外部磁場Ｇ2の下がり方が急激となっている。
【0241】
ユーザは磁場の掃引開始時刻と終了時刻、および量子相転移が起こる時刻とスケジューリング関数を指定すればよい。
図１２は、解析条件設定画面の一例を示す図である。例えば古典コンピュータ１００は、解析条件設定画面５０をモニタ２１に表示して、ユーザからの解析条件の設定入力を受け付ける。
図１２は、解析条件設定画面の一例を示す図である。例えば古典コンピュータ１００は、解析条件設定画面５０をモニタ２１に表示して、ユーザからの解析条件の設定入力を受け付ける。
【0242】
解析条件設定画面５０には、解析条件入力用のテキストボックス５１～５８と、実行指示用のボタン５９とが設けられている。テキストボックス５１は、掃引開始時刻ｔstartを設定する領域である。テキストボックス５２は、掃引終了時刻ｔendを設定する領域である。テキストボックス５３は、磁場のアニーリングスケジュールを変更する時刻ｔqptを設定する領域である。
【0243】
テキストボックス５４は、外部磁場Ｇ2の最大値Ｇ2,maxを設定する領域である。テキストボックス５４には、例えば初期値として「１．０」が設定されている。テキストボックス５５は、外部磁場Ｇ2の最小値Ｇ2,minを設定する領域である。テキストボックス５５には、例えば初期値として「０．０」が設定されている。テキストボックス５６は、磁場のアニーリングスケジュールを変更する時刻における外部磁場Ｇ2の値（Ｇ2,qpt）を設定する領域である。例えばテキストボックス５３に設定された値に応じて古典コンピュータ１００がＧ2,qptを計算し、計算結果をテキストボックス５６の初期値として設定してもよい。
【0244】
テキストボックス５７は、磁場のアニーリングスケジュールを変更する時刻以前に適用するアニーリングスケジュール関数の入力領域である。例えば使用可能な関数のリストがプルダウンメニュー５７ａに表示され、ユーザは、プルダウンメニュー５７ａから任意の関数を選択することができる。プルダウンメニュー５７ａ上で選択された関数が、テキストボックス５７に設定される。テキストボックス５８は、量子試行回数を設定する領域である。
【0245】
ユーザは、テキストボックス５１～５８に解析の条件を指定する情報を設定後、ボタン５９を押下することで、古典コンピュータ１００に対して、イジングモデルの基底状態を求める問題の解探索の実行を指示することができる。
【0246】
ここで大切なことは、量子相転移が起こるまではアニーリングスケジュールが重要であるが、一旦量子相転移が起こり、スピンが一斉に反転すると、後のアニーリングは高速に実行してよい点である。量子相転移が起こった時点で系の運命は決まっているからである。
【0247】
このため、重要なのはＧ20（つまり、ｔqpt）の値である。Ｇ20の値はスピンの磁気の強さＪに依存して決まる。スピンの磁気の強さＪを変えることで量子相転移が起こるＧ20の値は変わる。Ｇ20の値が０に近いほど計算には時間がかかる。Ｇ20の値が１に近いほど、高速に計算可能であるが、精度は良くない。そのため最適なＧ20の値が存在する。そのため、計算に都合の良いＧ20の値となるようにスピンの磁気の強さＪを設定することが求められる。
【0248】
適切なスピンの磁気の強さを求めるために、古典コンピュータ１００は、例えばスピンの強さＪを変えながら、一定回数、量子アニーリング型量子コンピュータ２００を制御してサンプリングを行い、その結果を記録する。サンプリングの結果としてスピン配置が得られる。古典コンピュータ１００は、サンプリング結果を用いてイジングモデルの全エネルギーが計算できる。そして、最低エネルギー状態が得られる頻度が得られる。その中で最低エネルギーが得られる頻度が最も高いＪの値を決めればよい。
【0249】
図１３は、スピンの磁気の強さ決定処理の手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。
［ステップＳ１０１］古典コンピュータ１００の処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最小値Ｊminの値を読み込む。
［ステップＳ１０１］古典コンピュータ１００の処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最小値Ｊminの値を読み込む。
【0250】
［ステップＳ１０２］処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最大値Ｊmaxの値を読み込む。
［ステップＳ１０３］処理部１２０は、記憶部１１０から横磁場強度を変えるサンプリング数Ｎesの値を読み込む。
［ステップＳ１０３］処理部１２０は、記憶部１１０から横磁場強度を変えるサンプリング数Ｎesの値を読み込む。
【0251】
［ステップＳ１０４］処理部１２０は、ループ回数を示す変数ｋが０からＮesになるまで、ｋの値をカウントアップしながらステップＳ１０５～Ｓ１０７の処理を繰り返す。
［ステップＳ１０５］処理部１２０は、スピンの磁気の強さを（３．２）式で求める。
［ステップＳ１０５］処理部１２０は、スピンの磁気の強さを（３．２）式で求める。
【0252】
【数75】
【0253】
［ステップＳ１０６］処理部１２０は、量子アニーリング型量子コンピュータ２００を用いたサンプリング処理を実行する。この際、処理部１２０は、サンプリング処理と同時に、最適なｔqptを決定することもできる。この処理の詳細は後述する（図１４～図１５参照）。
【0254】
［ステップＳ１０７］処理部１２０は、サンプリング処理によって得られた情報を、記憶部１１０に記録する。サンプリング処理で得られた情報は、スピン配列、最低エネルギー、および最低エネルギーに到達した頻度（確率）などである。
【0255】
［ステップＳ１０８］処理部１２０は、ループ回数ｋがＮesに達したら、処理をステップＳ１０９に進める。
［ステップＳ１０９］処理部１２０は、得られた頻度データから、一番低いエネルギー状態のスピン配置が最も高い頻度で得られたスピンの磁気の強さＪkの値を最適な横磁場強度として決定する。
［ステップＳ１０９］処理部１２０は、得られた頻度データから、一番低いエネルギー状態のスピン配置が最も高い頻度で得られたスピンの磁気の強さＪkの値を最適な横磁場強度として決定する。
【0256】
このようにして最適な横磁場強度が決定される。なお図１３に示すステップＳ１０６の処理では、量子アニーリング型量子コンピュータ２００を制御して、ｔqptも同時に決定することができる。例えば処理部１２０は、外部磁場（横磁場）の掃引開始時刻ｔstartと終了時刻ｔendを入力し、スケジュール関数ｓ（ｔ）を決めた上で、ｔqptを変えていく。そして処理部１２０は、最低エネルギー状態が得られる確率が最も高くなる磁場の掃引スケジュールにおけるｔqptを決める。
【0257】
以下図１４を参照して、ｔqptの決定を伴うサンプリング処理について詳細に説明する。
図１４は、ｔqptの決定を伴うサンプリング処理の手順の一例を示すフローチャートの前半である。以下、図１４に示す処理をステップ番号に沿って説明する。
図１４は、ｔqptの決定を伴うサンプリング処理の手順の一例を示すフローチャートの前半である。以下、図１４に示す処理をステップ番号に沿って説明する。
【0258】
［ステップＳ１２１］処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最小値Ｊminの値を読み込む。
［ステップＳ１２２］処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最大値Ｊmaxの値を読み込む。
［ステップＳ１２２］処理部１２０は、記憶部１１０からスピンの磁気の強さＪの最大値Ｊmaxの値を読み込む。
【0259】
［ステップＳ１２３］処理部１２０は、パラメータとして磁場掃引開始時刻ｔstartを設定する。例えば処理部１２０は、解析条件設定画面５０において入力された磁場掃引開始時刻ｔstartを設定する。
【0260】
［ステップＳ１２４］処理部１２０は、パラメータとして磁場掃引終了時刻ｔendを設定する。例えば処理部１２０は、解析条件設定画面５０において入力された磁場掃引終了時刻ｔendを設定する。
【0261】
［ステップＳ１２５］処理部１２０は、パラメータとして最大掃引回数ＮG2_loopを設定する。最大掃引回数ＮG2_loopは、例えば予めユーザによって指定されている。
［ステップＳ１２６］処理部１２０は、パラメータとして最大サンプリング回数Ｎsamplingを設定する。例えば処理部１２０は、解析条件設定画面５０において入力された量子試行回数の値を、最大サンプリング回数Ｎsamplingとして設定する。その後、処理部１２０は処理をステップＳ１３１（図１５参照）に進める。
［ステップＳ１２６］処理部１２０は、パラメータとして最大サンプリング回数Ｎsamplingを設定する。例えば処理部１２０は、解析条件設定画面５０において入力された量子試行回数の値を、最大サンプリング回数Ｎsamplingとして設定する。その後、処理部１２０は処理をステップＳ１３１（図１５参照）に進める。
【0262】
図１５は、ｔqptの決定を伴うサンプリング処理の手順の一例を示すフローチャートの後半である。以下、図１５に示す処理をステップ番号に沿って説明する。
［ステップＳ１３１］処理部１２０は、ループ回数を示す変数ｋが０からＮG2_loopになるまで、ｋの値をカウントアップしながらステップＳ１３２～Ｓ１３７の処理を繰り返す。
［ステップＳ１３１］処理部１２０は、ループ回数を示す変数ｋが０からＮG2_loopになるまで、ｋの値をカウントアップしながらステップＳ１３２～Ｓ１３７の処理を繰り返す。
【0263】
［ステップＳ１３２］処理部１２０は、量子相転移が起こる時刻ｔqptを（３．３）式によって設定する。
【0264】
【数76】
【0265】
［ステップＳ１３３］処理部１２０は、ループ回数を示す変数ｎが０からＮsamplingになるまで、ｎの値をカウントアップしながらステップＳ１３４～Ｓ１３６の処理を繰り返す。
【0266】
［ステップＳ１３４］処理部１２０は、初期状態の量子ビットに対してランダムなノイズを加える。
［ステップＳ１３５］処理部１２０は、量子アニーリング型量子コンピュータ２００に対して、ランダムなノイズが付加された量子ビットから、量子アニーリングによるイジングモデルの基底状態の探索処理を実行させる。処理部１２０は、探索処理の結果を、サンプリングデータとして量子アニーリング型量子コンピュータ２００から取得する。
［ステップＳ１３５］処理部１２０は、量子アニーリング型量子コンピュータ２００に対して、ランダムなノイズが付加された量子ビットから、量子アニーリングによるイジングモデルの基底状態の探索処理を実行させる。処理部１２０は、探索処理の結果を、サンプリングデータとして量子アニーリング型量子コンピュータ２００から取得する。
【0267】
［ステップＳ１３６］処理部１２０は、サンプリング結果に示されるスピン配列と、そのスピン配列におけるエネルギーとを記憶部１１０に記録する。
［ステップＳ１３７］処理部１２０は、ループ回数を示すｎの値がＮsamplingに達したら、処理をステップＳ１３８に進める。
［ステップＳ１３７］処理部１２０は、ループ回数を示すｎの値がＮsamplingに達したら、処理をステップＳ１３８に進める。
【0268】
［ステップＳ１３８］処理部１２０は、ループ回数を示すｋの値がＮG2_loopに達したら、処理をステップＳ１３９に進める。
［ステップＳ１３０］処理部１２０は、サンプリングデータの中からエネルギー最小値を求める。
［ステップＳ１３０］処理部１２０は、サンプリングデータの中からエネルギー最小値を求める。
【0269】
［ステップＳ１３１］処理部１２０は、エネルギー最小値の観測頻度が最も高いｔqptを、最低エネルギー状態が得られる確率が最も高くなる磁場の掃引スケジュールにおけるｔqptに決定する。
【0270】
このようにして、ｔqptが決定される。なお、最低エネルギー状態が得られる確率は、スピン単独の磁気の強さＪだけでなく、ノイズの与え方にも依存する。そのため最低エネルギー状態が得られる確率は、ユーザが明示的にノイズを与えた場合のノイズモデルにも依存し、量子アニーリング型量子コンピュータ２００が使用している量子プロセッサ（ＱＰＵ：Quantum Processing Unit）のノイズ特性にも依存する。従って、一般にはＱＰＵを変えれば、アニーリングの仕方はその都度変更することが求められる。
【0271】
図１６は、（１．１）式型の横磁場項と（２．１）式型の非線形横磁場を含んだ場合のＧ2に関する断熱ポテンシャルを示す図である。グラフ６１は、（１．１）型の線形の横磁場項を含んだ場合（非線形の外部磁場なし）の断熱ポテンシャルを示す。グラフ６２は（２．１）型の非線形の横磁場項を含んだ場合（非線形の外部磁場あり）の断熱ポテンシャルを示す。グラフ６１，６２の横軸はＧ2の値であり、縦軸はエネルギーである。グラフ６１，６２を比較すると、非線形磁場を含んだ場合の方が結果として得らえる量子状態の数が少ないことが分かる。
【0272】
ここで大切なことは量子相転移を起こすＧ2の値を制御することである。Ｇ20の値が１に近すぎればＧ2の値の変化量が大きい状況で量子相転移が起こる。この場合、最終的に得られるエネルギー値の分散は大きくなる。逆にＧ20の値が０に近すぎれば、今度は横磁場の効果が小さくなりすぎ、最終的に得られるエネルギー値の分散が大きくなる。このため、経験的にはＧ20の値が０．５≦Ｇ20≦０．８程度に収まるようにＧ20の値を調整するのが適切である。そのため処理部１２０は、横磁場の強さＪ＝Ｊiの値を変えながら、目標となるＧ20の値になるように調整する。
【0273】
＜虚時間発展に伴うエネルギー差の最適化＞
ここで、非線形の横磁場を導入した場合、なぜ、終状態として得られる量子状態の数が減るのかについて説明する。これは擬交差（avoided-crossing）の概念と、ゼーマン効果を理解することで簡単に説明がつく。
ここで、非線形の横磁場を導入した場合、なぜ、終状態として得られる量子状態の数が減るのかについて説明する。これは擬交差（avoided-crossing）の概念と、ゼーマン効果を理解することで簡単に説明がつく。
【0274】
図１７は、横磁場の強さＧ2に関するポテンシャルエネルギーの一例を示す図である。グラフ６３，６４の横軸はＧ2であり、縦軸はポテンシャルエネルギーである。グラフ６４は、グラフ６３の一部を拡大したものである。
【0275】
図１７の例では、簡単な例として３都市ＴＳＰが用いられている。３都市ＴＳＰの各都市の座標はそれぞれ（－１，０），（４，０），（０，２）である。このＴＳＰ問題をイジングモデルに直したとき、エネルギーの最小値が最小経路を与える。この問題の場合、最小エネルギーは５＋３√５≒１１．７０８である。この問題の解は１１．７０９である。
【0276】
図１７の例では、横磁場項を含めたハミルトニアンの行列形式に直して、対角化することでエネルギー固有値が求められている。具体的には基底関数（basis function）として、｜０００００００００＞から｜１１１１１１１１１＞で定義される２N個の基底関数を用いてハミルトニアン（１．１）式を直接対角化している。グラフ６３より、Ｇ2＝１からＧ2＝０まで連続的にエネルギーが変化し、Ｇ2＝０で問題の解である１１．７０８になっていることがわかる。
【0277】
グラフ６３，６４に示される曲線を、Ｇ2を断熱変数としたときの断熱ポテンシャルと呼ぶ。また、特定の断熱ポテンシャルのことを指してチャネルと呼ぶこともある。例えば、基底状態の断熱ポテンシャルを基底チャネルと呼ぶ場合もある。ここではエネルギーＥnに関連するチャネルをチャネル（ｎ）と呼ぶことにする。ｎ＝０の場合が基底状態である。
【0278】
図１７の例で用いた問題の場合、Ｇ2＝１において、第一励起状態にあっても基底状態に到達することができる。ここで、グラフ６４をよく見ると、断熱ポテンシャルが交差しているように見える部分がある。
【0279】
図１８は、断熱ポテンシャルが交差する部分を拡大した図である。グラフ６５は１段階拡大したものであり、グラフ６６は、グラフ６５よりもさらに大きく拡大したものである。
【0280】
グラフ６５において、Ｇ2＝０．９６７付近で断熱ポテンシャルが交差しているようにみえるが、拡大したグラフ６６を見ると分かるように実は交差していない。このように断熱ポテンシャルが交差しているように見えて、実は交差してない現象のことを疑交差と呼ぶ。疑交差は量子アニーリングでは励起状態への遷移を誘発する。つまり、疑交差は実質的に有限精度の計算では縮退しているようになってしまう。従って、疑交差が起こる領域ではＧ2に関するアニーリングスケジュールを遅くし、疑交差のエネルギー差ΔＥ1＝Ｅ1－Ｅ0よりも十分小さくなるようにＧ2を小さくすることが求められる。
【0281】
このエネルギー構造はハミルトニアン（１．１）式で決まってしまう。そのため、これを解消するために非線形の磁場項を導入するのである。
図１９は、非線形項として２次の横磁場を入れたハミルトニアン（２．１）式を対角化した例を示す図である。グラフ６７の横軸がＧ2であり、縦軸がエネルギーである。
図１９は、非線形項として２次の横磁場を入れたハミルトニアン（２．１）式を対角化した例を示す図である。グラフ６７の横軸がＧ2であり、縦軸がエネルギーである。
【0282】
グラフ６７を図１７の例と比較すると断熱ポテンシャルの構造がかなり違っていることがわかる。特にＧ2＝０．５付近では基底状態と第一励起状態のエネルギー差が大きくなっていることがわかる。Ｇ2＝０．５付近でエネルギー差が最大化するのはＧ3を（５．１）式のようにとったからである。
【0283】
【数77】
【0284】
（５．１）式ではＧ2＝０．５で非線形外部磁場が最大化するようになっているのだから、これは当然である。重要なことは外部磁場を導入することでエネルギー差を大きくすることが可能なことである。エネルギー差が大きくなれば、励起し辛くなる。これは量子アニーリング時に意図しない励起状態への励起が起こる確率を減らせることを意味する。ＴＳＰ３の場合は偽交差が起きても、基底状態に収束するような問題になっている。つまり、Ｇ2＝０で第一励起状態に属していたとしても、Ｅ＝１１．７０４の状態に収束することがわかる。しかし、４都市のＴＳＰにすると状況が一変する。
【0285】
図２０は、４都市ＴＳＰの場合のエネルギー準位と第ｎ励起状態のエネルギーＥnと基底状態のエネルギーＥ0との差を示す図である。Ｇ2＝０で基底状態だったものが、Ｇ2＝０．９付近で第一励起状態と疑交差することが分かる。グラフ６８が４都市ＴＳＰの場合のエネルギー準位を示している。グラフ６８の横軸がＧ2であり、縦軸がエネルギーである。グラフ６９が４都市ＴＳＰの場合の第ｎ励起状態のエネルギーＥnと基底状態のエネルギーＥ0との差を示している。グラフ６９の横軸がＧ2であり、縦軸がＥn－Ｅ0である。
【0286】
量子アニーリングにおいてＧ2＝０．９付近で虚数時間発展をしたときの収束解をＧ2＝０．９－δにおける初期条件にすると、δの値が十分小さくないと、初期状態のエネルギーがこの二つの状態よりも高くなってしまう。この状態で虚数時間発展をすると、Ｅ＝７．４１４の量子状態にたどり着くか、それとも、Ｅ＝７．８１２の状態にたどり着くかはランダムになってしまう。
【0287】
この場合、Ｅ＝７．４１４の量子状態と、Ｅ＝７．８１２の量子状態のどちらにたどり着くかを決める分水嶺はグラフ６９から分かる。疑交差が生じるＧ2の領域で、ΔＥ1とΔＥ2が２つに分岐している。一方はＥ＝７．４１４に到達する。もう一方は、Ｅ＝７．８１２に到達する。グラフ６９は虚数時間発展により基底状態に到達することを保証するエネルギー差に関する条件を与える図になっている。あるＧ2に対して、虚数時間発展を行い、得られた終状態をＧ2－δにおける初期状態とする。このとき、この初期状態のエネルギーはわずかに上がることになるが、その値が、ΔＥ1が７．４１４に収束する領域内に入っていれば、計算は基底状態に収束する。つまり、Ｇ2を変える量子アニーリングスケジュールは、Ｇ2－δにおける初期状態のエネルギーが一定値以下になるように、エネルギー変化量から決めるべきものである。
【0288】
従って処理部１２０は、Ｇ2を変えるときに、図１３に示したフローチャートのようにエネルギーを計算しながら、指定した閾値よりも小さくなるようにＧ2を決める。なお、グラフ６９は点線が線形の横磁場の結果であり、実線が非線形磁場を含んだ結果である。非線形磁場を含んだ場合、エネルギー差は線形磁場より常に大きくなっているため、収束領域が増す効果がある。
【0289】
従って、Ｇ2を大きく変えすぎるとエネルギーが大きく変わってしまうため、容易にエネルギーが高いチャネルへ遷移してしまう。エネルギーが高いチャネルへ遷移してしまうと、Ｇ2を小さくしていったとき、量子相転移が起こったときに、最低エネルギー状態とは違うスピン配置になるようにスピンが反転してしまう。このため、Ｇ2の値の弱め方は予め閾値ε1を決め、Ｇ2－δにおける終状態のエネルギーＥ（Ｇ2－δ）と、Ｇ2における終状態のエネルギーＥ（Ｇ2）との差ΔＥ＝｜Ｅ（Ｇ2－δ）－Ｅ（Ｇ2）｜≦ε1を満たすようにして計算をしていくことが求められる。
【0290】
図２１は、Ｇ2の弱め方の手順の一例を示すフローチャートである。以下、図２１に示す処理をステップ番号に沿って説明する。
［ステップＳ２０１］処理部１２０は、ループ回数を示す変数ｎが１からＮloopになるまで、ｎの値をカウントアップしながらステップＳ２０２～Ｓ２１６の処理を繰り返す。Ｎloopは予め与えられた自然数である。
［ステップＳ２０１］処理部１２０は、ループ回数を示す変数ｎが１からＮloopになるまで、ｎの値をカウントアップしながらステップＳ２０２～Ｓ２１６の処理を繰り返す。Ｎloopは予め与えられた自然数である。
【0291】
［ステップＳ２０２］処理部１２０は、Ｇ2,n（ｎ回目の量子試行におけるＧ2）の値を読み込む。
［ステップＳ２０３］処理部１２０は、δの値の初期値を読み込む。
［ステップＳ２０３］処理部１２０は、δの値の初期値を読み込む。
【0292】
［ステップＳ２０４］処理部１２０は、Ｇ'2,n＝Ｇ2,n－δとして、Ｇ'2,nを初期化する。
［ステップＳ２０５］処理部１２０は、エネルギー差の閾値ε1を読み込む。閾値ε1は、予め設定された実数である。
［ステップＳ２０５］処理部１２０は、エネルギー差の閾値ε1を読み込む。閾値ε1は、予め設定された実数である。
【0293】
［ステップＳ２０６］処理部１２０は、ループ回数を示す変数ｋが１からＮbisectionになるまで、ｋの値をカウントアップしながらステップＳ２０７～Ｓ２１５の処理を繰り返す。Ｎbisectionは予め与えられた自然数である。
【0294】
［ステップＳ２０７］処理部１２０は、Ｇ2,n-1の値での終状態Ｅn-1,finalの｜Ψn-1,final＞を｜Ψn,initial＞とする。
［ステップＳ２０８］処理部１２０は、｜Ψn,initial＞を用いて初期エネルギーＥn,initial＝｜Ψn,initial｜Ｈ｜Ψn,initial＞（Ｈは＾付き）を計算する。
［ステップＳ２０８］処理部１２０は、｜Ψn,initial＞を用いて初期エネルギーＥn,initial＝｜Ψn,initial｜Ｈ｜Ψn,initial＞（Ｈは＾付き）を計算する。
【0295】
［ステップＳ２０９］処理部１２０は、ΔＥ＝｜Ｅn-1,final－Ｅn,initial｜を計算する。
［ステップＳ２１０］処理部１２０は、虚数時間を発展させ、Ｇ'2,nにおける最低エネルギー状態を求める。
［ステップＳ２１０］処理部１２０は、虚数時間を発展させ、Ｇ'2,nにおける最低エネルギー状態を求める。
【0296】
［ステップＳ２１１］処理部１２０は、Ｇ’2,nの値での終状態Ｅn,finalおよび｜Ψn-1,final＞を保存する。
［ステップＳ２１２］処理部１２０は、ΔＥ’＝｜＜Ψn,final｜Hn｜Ψn,final＞－Ｅn,initial｜（Ｈは＾付き）を計算する。
［ステップＳ２１２］処理部１２０は、ΔＥ’＝｜＜Ψn,final｜Hn｜Ψn,final＞－Ｅn,initial｜（Ｈは＾付き）を計算する。
【0297】
［ステップＳ２１３］処理部１２０は、ΔＥ≦ε1かつΔＥ’≦ε1の条件が満たされるか否かを判断する。処理部１２０は、条件が満たされる場合、処理をステップＳ２１６に進める。また処理部１２０は、条件が満たされない場合、処理をステップＳ２１４に進める。
【0298】
［ステップＳ２１４］処理部１２０は、δの値をδ／２に更新する（δ＝δ／２）。
［ステップＳ２１５］処理部１２０は、Ｇ’2,nをＧ2,n－δにする（Ｇ’2,n＝Ｇ2,n－δ）。
［ステップＳ２１５］処理部１２０は、Ｇ’2,nをＧ2,n－δにする（Ｇ’2,n＝Ｇ2,n－δ）。
【0299】
［ステップＳ２１６］処理部１２０は、変数ｋの値がＮbisectionに達したら処理をステップＳ２１７に進める。
［ステップＳ２１７］処理部１２０は、変数ｎの値がＮloopに達したら処理を終了する。
［ステップＳ２１７］処理部１２０は、変数ｎの値がＮloopに達したら処理を終了する。
【0300】
このようにして、励起チャネルに量子状態が遷移しないようにＧ2を小さくすることができる。
なお図９に示したような結果から、非線形横磁場を加えることで、基底状態含め、低エネルギー状態が得られる確率が増えたことは分かる。ただし、基底状態のみが得られるわけではない。その原因は、横磁場項の次数に関連する。
なお図９に示したような結果から、非線形横磁場を加えることで、基底状態含め、低エネルギー状態が得られる確率が増えたことは分かる。ただし、基底状態のみが得られるわけではない。その原因は、横磁場項の次数に関連する。
【0301】
この原理を説明するために、４スピン系の例題を出して説明する。今、簡単のため、４スピン系の｜００００＞という状態を考える。この量子状態に虚数時間の時間発展演算子ｅｘｐ｛－ＴΔｔ｝（Ｔは＾付き）を乗じると、（５．２）式となる。
【0302】
【数78】
【0303】
横磁場演算子を乗じることで、１スピンが反転した｜１０００＞のような状態が得られる。２回目を乗じると（５．３）式となる。
【0304】
【数79】
【0305】
２粒子励起が出てくるように思えるが、実は出ない。支配項は相変わらず、｜００００＞であるし、Δｔに比例する１粒子励起の項の係数が演算回数分だけ増幅されているだけである。つまり、「２粒子励起は実質上無視される」こととなる。
【0306】
次に、作用素をＴ＋Ｔ2（Ｔは共に＾付き）にしたらどうなるかを考察する。鈴木－トロッター分解は１次の範囲で（５．４）式となる。
【0307】
【数80】
【0308】
従って（５．５）式が得られる。
【0309】
【数81】
【0310】
そのため（５．６）式となる。
【0311】
【数82】
【0312】
ここからわかる通り、高次の磁場を導入することで２粒子励起以上の項をΔｔの項に含めることができる。そのため、非線形磁場を加えたときの遷移確率は、非線形磁場が摂動的であれば、（５．７）式で与えられる。
【0313】
【数83】
【0314】
ここにρは終状態の状態密度である。これはフェルミの黄金律（Fermi’s golden rule）と呼ばれる。つまり、外部から非線形磁場を加えることで２スピン以上が関わった断熱チャネル間で結合が起こる。結合が起きれば、結合強度に応じて縮退が解ける。そのため、基底状態と結合したチャネルが観測されなくなっているのである。これは励起状態でも同じである。励起状態同士で結合した場合、結合したチャネルの中でエネルギーが一番低いモードが観測されることになる。
【0315】
図２２は、４都市ＴＳＰのハミルトニアンを数値的に対角化したときの固有ベクトルの一例を示す図である。図２２に示すグラフ７１～７４は、横軸が固有状態それぞれを区別するインデックスであり、縦軸が各固有状態に対応する固有ベクトルの要素値である。グラフ７１は、Ｇ2＝１．０の場合である。横磁場項の基底状態のため、全ての状態が等しい重さで入っている。グラフ７２は、Ｇ2＝０．９９の場合である。縦磁場項（求解ハミルトニアン）の情報が入るため、ここで基底に重みが加わる。そして、Ｇ2の値を小さくしていくと、特定の基底の振幅が増幅していく（グラフ７３，７４）。
【0316】
グラフ７４に示すように最終的にはこの場合２つの量子状態のみが残るが、これは両方ともＥ＝７．４１４の状態である。基底状態が２つだけ求まっているように見えるのは、最初の２０４８次元分だけを取り出しているからである。
【0317】
図２３は、固有ベクトルのインデックスと固有状態との対応関係の一例を示す図である。図２３に示す対応表８１は、図２２のグラフ７４において、固有ベクトルの振幅が大きいインデックスを幾つか抜粋し、対応するスピン状態、エネルギー固有値をまとめたものである。スピン状態を見れば一目瞭然であるが、ハミング距離（Hamming Distance）が１のものは有限の振幅を持っているが、２以上のものは急激に固有ベクトルの振幅が小さくなっていることが分かる。つまり、これは横磁場を１次にしているからである。基底状態に対してハミング距離２以上のモードは実質上無視されることになる。
【0318】
この事実は次のことを暗示する。２次の横磁場を入れた場合、ある量子状態｜ｉ＞と結合する（遷移可能）であるのは、高々ハミング距離が２以下、つまり、スピン２個が反転したモードまでである。一般に、基底状態、第一励起状態のスピン状態のハミング距離は大きく違うため、２次の横磁場を入れれば、それぞれの状態からハミング距離２以下のモードが結合する。そして、一種のグループ化が起こる。結果として、観測される量子状態はグループ化が起こったため減る。
【0319】
非線形磁場を導入することで複数の系列の断熱ポテンシャルのエネルギー差が大きくなるが、互いに影響を及ぼさないチャネル同士のエネルギー差は非線形横磁場項の影響を受けない。それでも、部分的には影響を受けるため、最終的に観測される量子状態が減るのである。ただし、２次の非線形の横磁場は高々２スピンが関与した素過程しか記述しないため、必ずしもすべてのチャネルが結合に関与するわけではない。このため、基底状態のみが観測される状態にはできないのである。
【0320】
そして、この事実は次の事実も暗示する。一般にすべての状態を結合させるための横磁場の次数はＮである。横磁場の次数は１次項からＮ次まですべて加えることが求められる。あるいは、ｅｘｐ（－βＴ）（Ｔは＾付き）のように、級数展開がＮ次まで定義できるような非線形横磁場を加えるべきである。
【0321】
ただし、２次の横磁場を用いても、部分的には縮退が解けるのだから、基底状態が得られる確率が向上するため、非常に簡単な非線形磁場を用いても量子アニーリングの効率は向上する。
【0322】
また、２次の横磁場を用いた場合に、特定の構造を持つハミルトニアン系において、基底状態が得られる確率が劇的に向上する場合もありえる。これはハミルトニアンの構造に依存するが、２次の非線形横磁場で結合するモードが非常に多くなれば、それだけ効率は上がることになる。
【0323】
以上煩雑な説明をしたが、現象の本質は簡単である。初等的な量子力学の教科書に書いてある「ゼーマン効果」なのである。あるハミルトニアン系で固有値に縮退が起こっていた。磁場を加えてハミルトニアン行列の非対角項を導入することで縮退が解けるということである。縮退が解けるのは、非線形外部磁場によりハミルトニアン行列の非対角項が０でない者同士の組に限られる。そのため、互いに影響を及ぼさないグループが存在する。影響を及ぼしあう、つまり、結合するチャネル同士はグループ化される。結果として、終状態として観測される状態数は減少する。結果として基底状態が得られる確率が相対的に上がるのである。もし、縮退を完全に解きたいのであれば、Ｎスピンを反転させることが求められる。そのため、（５．８）式のような横磁場項を用いることとなる。
【0324】
【数84】
【0325】
そして、Ｇ2＝０とＧ2＝１における境界条件を満たすためには、Ｇ2＝０とＧ2＝１においてｇ（Ｇ2）＝０を満たす関数ｇ（Ｇ2）を用いて、ハミルトニアンを（５．９）式とすることとなる。
【0326】
【数85】
【0327】
以上のようにして、ノイズを付加した初期状態から、非線形横磁場を追加して、実数時間発展と虚数時間発展とを組合せて量子アニーリングを行い、アニーリングスケジュールを適切なＧ2を境に変更し、計算精度を落とさずに処理の効率化を図ることができる。
【0328】
＜量子相転移の発生有無の検知＞
図１２に示す様な解析条件設定画面５０で量子相転移が発生する時点をｔqptまたはG2,qptで指定可能であるが、アニーリングの過程で量子相転移の発生の有無を検知し、量子相転移が発生した後、アニーリングのスケジュールを変更することもできる。例えば古典コンピュータ１００においてＳＡを行えばアニーリングの過程でのエネルギーを計算できる。そこで虚数時間発展の過程でのエネルギー遷移に基づいて、量子相転移の発生を検知することができる。
図１２に示す様な解析条件設定画面５０で量子相転移が発生する時点をｔqptまたはG2,qptで指定可能であるが、アニーリングの過程で量子相転移の発生の有無を検知し、量子相転移が発生した後、アニーリングのスケジュールを変更することもできる。例えば古典コンピュータ１００においてＳＡを行えばアニーリングの過程でのエネルギーを計算できる。そこで虚数時間発展の過程でのエネルギー遷移に基づいて、量子相転移の発生を検知することができる。
【0329】
図２４は、量子相転移の検出を伴うアニーリング処理の手順の一例を示すフローチャートである。以下、図２４に示す処理をステップ番号に沿って説明する。
［ステップＳ３０１］処理部１２０は、スピンの磁気の強さＪの値を読み込む。
［ステップＳ３０１］処理部１２０は、スピンの磁気の強さＪの値を読み込む。
【0330】
［ステップＳ３０２］処理部１２０は、Ｇ2,min＝０，Ｇ2,max＝１を設定する。
［ステップＳ３０３］処理部１２０は、ループ回数を示す変数ｎの値を１ずつ減算しながら、ＮG2_loopから０まで、ステップＳ３０４～Ｓ３０７の処理を繰り返す。
［ステップＳ３０３］処理部１２０は、ループ回数を示す変数ｎの値を１ずつ減算しながら、ＮG2_loopから０まで、ステップＳ３０４～Ｓ３０７の処理を繰り返す。
【0331】
［ステップＳ３０４］処理部１２０は、横磁場の強度をべき乗のアニーリングスケジュール「Ｇ2,k＝（ｎ／ＮG2_loop）-c」を用いてＧ2,kを決定する。
［ステップＳ３０５］処理部１２０は、虚数時間発展の処理を実行する。
［ステップＳ３０５］処理部１２０は、虚数時間発展の処理を実行する。
【0332】
［ステップＳ３０６］処理部１２０は、虚数時間発展中にエネルギーの２回微分「ｄ2Ｅn／ｄＧ2,n
2」を計算する。
［ステップＳ３０７］処理部１２０は、「（ｄ2Ｅk／ｄＧ2,k 2）（ｄ2Ｅk-1／ｄＧ2,k-1 2）＜０」の条件が満たされるか否かを判断する。処理部１２０は、条件が満たされる場合、処理をステップＳ３０９に進める。また処理部１２０は、条件が満たされない場合、処理をステップＳ３０８に進める。
［ステップＳ３０７］処理部１２０は、「（ｄ2Ｅk／ｄＧ2,k 2）（ｄ2Ｅk-1／ｄＧ2,k-1 2）＜０」の条件が満たされるか否かを判断する。処理部１２０は、条件が満たされる場合、処理をステップＳ３０９に進める。また処理部１２０は、条件が満たされない場合、処理をステップＳ３０８に進める。
【0333】
［ステップＳ３０８］処理部１２０は、変数ｎの値が０に達した場合、処理をステップＳ３０９に進める。
［ステップＳ３０９］処理部１２０は、量子相転移が起こったと判断し、この時点のＧ2,kの値を保存する。
［ステップＳ３０９］処理部１２０は、量子相転移が起こったと判断し、この時点のＧ2,kの値を保存する。
【0334】
［ステップＳ３１０］処理部１２０は、Ｇ2,min＝０，Ｇ2,max＝Ｇ2,k，Ｎ’G2_loopを設定する。Ｎ’G2_loopは、例えば現在のｎよりも小さい値とする（Ｎ’G2_loop＜ｎ）。
［ステップＳ３１１］処理部１２０は、、ループ回数を示す変数ｎ’の値を１ずつカウントアップしながら、０からＮ’G2_loopまで、ステップＳ３１２～Ｓ３１３の処理を繰り返す。
［ステップＳ３１１］処理部１２０は、、ループ回数を示す変数ｎ’の値を１ずつカウントアップしながら、０からＮ’G2_loopまで、ステップＳ３１２～Ｓ３１３の処理を繰り返す。
【0335】
［ステップＳ３１２］処理部１２０は、横磁場の強度を「Ｇ2＝Ｇ2,k－ｎ’（Ｇ2,k／Ｎ’G2_loop－１）」に設定する。
［ステップＳ３１３］処理部１２０は、虚数時間発展の処理を実行する。
［ステップＳ３１３］処理部１２０は、虚数時間発展の処理を実行する。
【0336】
［ステップＳ３１４］処理部１２０は、変数ｎ’の値がＮ’G2_loopに達した場合、処理を終了する。
このように虚数時間発展のなかで量子相転移の発生を検知して、アニーリングスケジュールを変更することが可能となる。アニーリングスケジュールの変更後は、変更前に予定していた残りのサンプリング回数よりも少ないサンプリング回数で、イジングモデルの終状態が得られる。これにより、処理の効率化を図ることができる。
このように虚数時間発展のなかで量子相転移の発生を検知して、アニーリングスケジュールを変更することが可能となる。アニーリングスケジュールの変更後は、変更前に予定していた残りのサンプリング回数よりも少ないサンプリング回数で、イジングモデルの終状態が得られる。これにより、処理の効率化を図ることができる。
【0337】
以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。
【符号の説明】
【0338】
１ 量子アニーリング型量子コンピュータ
１０ 情報処理装置
１１ 記憶部
１２ 処理部
１０ 情報処理装置
１１ 記憶部
１２ 処理部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】