特許 7518379 最適化装置、最適化プログラム、および最適化方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-07-09

(45)【発行日】2024-07-18

(54)【発明の名称】最適化装置、最適化プログラム、および最適化方法

(51)【国際特許分類】

   G06N 10/60 20220101AFI20240710BHJP

   G06N 99/00 20190101ALI20240710BHJP

【ＦＩ】

G06N10/60

G06N99/00 180

【請求項の数】 9

(21)【出願番号】P 2020170876

(22)【出願日】2020-10-09

(65)【公開番号】P2022062760

(43)【公開日】2022-04-21

【審査請求日】2023-07-07

(73)【特許権者】

【識別番号】000005223

【氏名又は名称】富士通株式会社

(74)【代理人】

【識別番号】110002918

【氏名又は名称】弁理士法人扶桑国際特許事務所

(72)【発明者】

【氏名】櫛部 大介

【審査官】多賀 実

(56)【参考文献】

【文献】特開２０１８－０１０４７４（ＪＰ，Ａ）

【文献】特開２０１９－０７９２２５（ＪＰ，Ａ）

【文献】特開２０２０－００９３０１（ＪＰ，Ａ）

【文献】米国特許出願公開第２０１４／０２９７２４７（ＵＳ，Ａ１）

【文献】米国特許出願公開第２０２０／００９００７１（ＵＳ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｎ ３／００－９９／００

(57)【特許請求の範囲】

【請求項1】

求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化のシミュレーションをすることで前記イジングモデルの基底状態を求解する際に、前記シミュレーションで使用する係数の一部にノイズに相当する値を付加し、前記磁場の強度を前記シミュレーション上の時間の進行に伴って低減させていく実数時間発展の第１処理と、虚数時間発展法に基づき、前記イジングモデルのエネルギーを低減させる第２処理とを実行する処理部、
を有する最適化装置。

【請求項2】

前記処理部は、前記イジングモデルの初期状態に対して前記ノイズに相当する値を付加し、前記ノイズに相当する値が付加された前記イジングモデルの状態を前記シミュレーション開始時の状態とする、
請求項１記載の最適化装置。

【請求項3】

前記処理部は、前記イジングモデルに含まれる複数のスピンのうちの一部のスピンの波動関数に含まれる係数に前記ノイズに相当する値を付加する、
請求項２記載の最適化装置。

【請求項4】

前記処理部は、前記ノイズに相当する値を付加する前記一部のスピンを、前記複数のスピンの中からランダムに選択する、
請求項３記載の最適化装置。

【請求項5】

前記処理部は、前記イジングモデルに対する前記ノイズに相当する値の付加を複数回実行することで、異なる前記ノイズが付加された前記イジングモデルの複数の状態を生成し、前記イジングモデルの前記複数の状態それぞれを前記シミュレーション開始時の状態として前記第１処理と前記第２処理とを実行する、
請求項２ないし４のいずれかに記載の最適化装置。

【請求項6】

前記処理部は、前記イジングモデルに含まれるスピンにおける前記磁場の強さを示す係数に、前記ノイズに相当する値を付加する、
請求項１記載の最適化装置。

【請求項7】

前記処理部は、前記第１処理と前記第２処理とを繰り返し実行することで得られた、量子力学に基づく前記イジングモデルのスピンの状態を、古典力学に基づいて取り得る値のいずれかに決定し、各スピンが決定された値を有する古典イジングモデルの状態を初期状態として、古典力学に基づく解法によって前記イジングモデルの基底状態を求める第３処理を実行する、
請求項１ないし６のいずれかに記載の最適化装置。

【請求項8】

コンピュータに、
求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化のシミュレーションをすることで前記イジングモデルの基底状態を求解する際に、前記シミュレーションで使用する係数の一部にノイズに相当する値を付加し、前記磁場の強度を前記シミュレーション上の時間の進行に伴って低減させていく実数時間発展の第１処理と、虚数時間発展法に基づき、前記イジングモデルのエネルギーを低減させる第２処理と、
を実行させる最適化プログラム。

【請求項9】

コンピュータが、
求解対象の問題を表すイジングモデルに印加される磁場を低減させたときの前記イジングモデルの状態変化のシミュレーションをすることで前記イジングモデルの基底状態を求解する際に、前記シミュレーションで使用する係数の一部にノイズに相当する値を付加し、前記磁場の強度を前記シミュレーション上の時間の進行に伴って低減させていく実数時間発展の第１処理と、虚数時間発展法に基づき、前記イジングモデルのエネルギーを低減させる第２処理と、
を実行する最適化方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、最適化装置、最適化プログラム、および最適化方法に関する。

【背景技術】

【0002】

自然科学や社会科学において頻出する問題として、評価関数（エネルギー関数とも呼ばれる）の最小値（または最小値を与える評価関数の状態変数の値の組合せ）を求める最小値求解問題（または組合せ最適化問題と呼ばれる）がある。近年、このような問題を磁性体のスピンの振る舞いを表すモデルであるイジングモデルで定式化する動きが加速している。

【0003】

よく使用される解法は、マルコフ連鎖モンテカルロ（ＭＣＭＣ：Markov Chain Monte Carlo）法を用い確率過程を導入し、ボルツマン分布など特定の分布を定義したうえで温度を導入し、温度をシミュレーション的に下げる方法である。この解法はシミュレーテッドアニーリング（ＳＡ：Simulated Annealing）法と呼ばれる。しかし、ＳＡ法で厳密解に到達するためには温度スケジュールを時間に関して対数の逆数でとることとなる。すなわち、高速に温度の冷却を行うと、必ずしも解が最適解に達しない。

【0004】

このような状況下で、近年量子力学に基づく計算技術および計算機の開発の動きが活発化している。この動きの技術的な基盤はイジング型量子コンピュータの実現である。イジング型の量子コンピュータは量子アニーリング（ＱＡ：Quantum Annealing）法が理論的な基盤になっている。量子コンピュータは普段使用している電子回路をベースにした古典的なコンピュータでは現実的な計算時間内で解くことが難しい問題を解くことができると期待されている。現段階ではすべての問題が量子コンピュータで高速化できるわけではないが、一部の問題が量子コンピュータで高速に解けると期待されている。代表的な問題がＲＳＡ暗号など暗号通信の基礎になっている素因数分解の問題である。このように、イジング型ではあるが量子コンピュータが実現されたため、量子コンピュータの応用範囲を探索する動きが加速している。

【先行技術文献】

【特許文献】

【0005】

【文献】特開２０１８－０６７２００号公報

【文献】特開２０２０－６４５３５号公報

【文献】特開２０２０－９３０１号公報

【非特許文献】

【0006】

【文献】Sanroku Tsukamoto, Motomu Takatsu, Satoshi Matsubara and Hirotaka Tamura, “An Accelerator Architecture for Combinatorial Optimization Problems”, FUJITSU Sci. Tech. J. Vol.53, No.5, pp.8-13, September 2017

【文献】S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi,“Optimization by Simulated Annealing”, Science, Vol.220, No.4598, pp.671-680, 13, May, 1983

【文献】Koji Hukushima and Koji Nemoto, “Exchange Monte Carlo Method and Application to Spin Glass Simulations”, Journal of the Physical Society of Japan, Volume 65, No.6, pp.1604-1608, June, 1996

【文献】E. D. Dahl, “Programming with D-Wave: Map Coloring Problem”, D-Wave Systems whitepaper, November 2013

【文献】Tadashi Kadowaki and Hidetoshi Nishimori, “Quantum annealing in the transverse Ising model”, Phys. Rev. E, 58, No.5, 5355, November 1998

【文献】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

【文献】T. Albash and D. A. Lidar, “Demonstration of a scaling advantage for a quantum annealer over simulated annealing”, arXiv:1705.07452v3, 6 August 2018

【文献】J. R. Abo-Shaeer, C. Raman, J. M. Vogels, W. Ketterle, “Observation of Vortex Lattices in Bose-Einstein Condensates”, Science, Vol. 292, Issue 5516, pp.476-p479, 20 April 2001

【文献】Masaki Mutou, Toru Morishita, Daisuke Kushibe, Shinichi Watanabe, and Michio Matsuzawa, “Numerical simulation of quantized vortices in 2-Dimensional Bose-Einstein Condensates”, The Fifth Asian International Seminar on Atomic and Molecular Physics , 202-203, October 2002

【文献】D. L. Feder, M.S. Pindzola, L. A. Collins, B. I. Schneider and C. W. Clark, “Dark-soliton states of Bose-Einstein condensates in anisotropic traps”, Phys. Rev. A, 62, 053606(2000), 26 Apr 2000

【文献】Naomichi Hatano and Masuo Suzuki, “Finding Exponential Product Formulas of Higher Orders”, Quantum Annealing and Other Optimization Methods, pp.37-68, 16 November 2005

【文献】中原幹夫、坂東将光、田中宗、“断熱量子コンピューティングによる巡回セールスマン問題の解法”近畿大学理工学総合研究所報告29、1-9, 28 February 2017

【発明の概要】

【発明が解決しようとする課題】

【0007】

Ｎ個（Ｎは１以上の整数）の２値変数で定義された最小値求解問題において、量子力学に原理を置く計算を行い系の基底状態（最適値）および基底状態の波動関数（最適値における状態変数の組み合わせ）を求める問題を考える。ここで系とは、周囲の環境とは切り離して考えた部分空間を指す。この問題をイジングモデルに変換すれば、ＱＡ法におけるハミルトニアンを用いて解くことができる。最小値求解問題を求解する装置は、最適化装置と呼ばれる。

【0008】

最適化装置は、ＱＡ法を用いた解探索では、初期状態としてＮ個のスピンに横磁場をかけた状態にしておく。また最適化装置は、初期状態として実験的に準備可能な量子状態を用意する。最適化装置は、その状態から、時間経過と共に横磁場を弱めていく。すると、最終的に横磁場が完全になくなり、求解対象のハミルトニアンの基底状態、すなわちエネルギー最小値が求まる。

【0009】

しかし、ＱＡ法を用いた解探索では、１次の量子相転移と呼ばれる現象の発生を抑制するため、横磁場の弱め方には強い制限がかかっている。すなわち最適化装置は、横磁場の弱め方を遅くすることで１次の量子相転移の発生を抑止することができる。しかしながら、１次の量子相転移の発生を抑止できるように横磁場の弱め方を遅くすると、最小値求解問題の解探索時間が長期化する。

【0010】

１つの側面では、本件は、イジングモデルの基底状態を効率的に求めることを目的とする。

【課題を解決するための手段】

【0011】

１つの案では、以下の処理を実行する処理部を有する最適化装置が提供される。
処理部は、求解対象の問題を表すイジングモデルに印加される磁場を低減させたときのイジングモデルの状態変化のシミュレーションをすることでイジングモデルの基底状態を求解する際に、シミュレーションで使用する係数の一部にノイズに相当する値を付加する。そして処理部は、磁場の強度をシミュレーション上の時間の進行に伴って低減させていく実数時間発展の第１処理と、虚数時間発展法に基づき、イジングモデルのエネルギーを低減させる第２処理とを実行する。

【発明の効果】

【0012】

１態様によれば、イジングモデルの基底状態を効率的に求めることができる。

【図面の簡単な説明】

【0013】

【図1】第１の実施の形態に係る最適化方法の一例を示す図である。

【図2】最適化装置のハードウェアの一例を示す図である。

【図3】虚数時間発展法を用いた熱的散逸の模式図である。

【図4】４都市のＴＳＰにおける都市の配置例を示す図である。

【図5】Ｇ₂の値の数に応じた計算結果の違いを示す図である。

【図6】横磁場の強さに応じたエネルギー値を模式的に示す図である。

【図7】ノイズ付加を伴うＱＴＳＡ計算を模式的に示す図である。

【図8】最適化装置の機能を示すブロック図である。

【図9】解探索手順の一例を示すフローチャートである。

【図10】ハミルトニアン情報読み込み処理の手順の一例を示すフローチャートである。

【図11】スピン情報初期化処理の手順の一例を示すフローチャートである。

【図12】量子試行計算処理の手順の一例を示すフローチャートである。

【図13】波動関数へのノイズ付加処理の手順の一例を示すフローチャートである。

【図14】虚数時間発展を用いたＱＡ処理の一例を示すフローチャートである。

【図15】対応する古典的状態を求める処理の手順の一例を示すフローチャートである。

【図16】４都市ＴＳＰにおける２００回の量子試行結果の一例を示す図である。

【図17】スピンの波動関数の時間変化の一例を示す図である。

【図18】８都市のＴＳＰの一例を示す図である。

【図19】８都市のＴＳＰの量子試行結果の一例を示す図である。

【図20】８都市ＴＳＰにおける波動関数の時間発展の一例を示す図である。

【図21】８都市ＴＳＰにおける量子試行と古典イジングモデルのＭＣＭＣとの計算結果の比較例を示す図である。

【図22】試行回数に応じたエネルギーの遷移例を示す図である。

【図23】ulysses16におけるＱＴＳＡと古典ＭＣＭＣで得られた解の比較例を示す図である。

【図24】解の改善速度の量子論と古典論との比較例を示す図である。

【図25】典型的な量子試行の時間発展の様子を示す図である。

【図26】量子論におけるエネルギー散逸の様子を示す図である。

【図27】ulysses22におけるランキング情報の一例を示す図である。

【図28】ulysses22における解の更新の様子を示す図である。

【図29】エネルギー分解能ごとの量子試行結果の一例を示す図である。

【図30】横磁場にノイズを入れた場合のＱＴＳＡによる量子アニーリング結果の一例を示す図である。

【発明を実施するための形態】

【0014】

以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
第１の実施の形態は、イジングモデルの基底状態の求解において熱的散逸の効果を取り入れるため、量子力学的な時間発展演算子における時間を複素数に拡張し、虚数時間発展法に基づき、熱的散逸を実現する最適化方法である。熱的散逸機構の導入により、断熱遷移条件を厳密に守らずとも基底状態へ緩和可能である。なお虚数時間発展法を用いて熱的散逸機構を導入しても、基底状態が求まらない場合があり得る。そのため第１の実施の形態に係る最適化方法では、波動関数の対象性を破壊するノイズをイジングモデルに付加することにより、基底状態に到達する可能性を向上させている。

【0015】

図１は、第１の実施の形態に係る最適化方法の一例を示す図である。図１には、イジングモデル１のエネルギーの最適化方法を実施する最適化装置１０が示されている。最適化装置１０は、例えば最適化プログラムを実行することにより、最適化方法を実施することができる。

【0016】

最適化装置１０は、記憶部１１と処理部１２とを有する。記憶部１１は、例えば最適化装置１０が有するメモリ、またはストレージ装置である。処理部１２は、例えば最適化装置１０が有するプロセッサ、または演算回路である。

【0017】

記憶部１１は、求解対象となる問題を示す求解問題情報１１ａを記憶する。求解対象となる問題は、例えばイジングモデル１に変換可能な最小値求解問題である。求解問題情報１１ａには、例えば求解対象の問題を示すイジングモデル１の情報が設定される。イジングモデル１は２値Ｎ変数（２値の変数がＮ個あること）で定義され、イジングモデル１の最終的な基底状態が、求解対象の問題の最適解を表す。

【0018】

処理部１２は、量子力学に基づく変分原理を用いてイジングモデル１の基底状態を求解する。すなわち処理部１２は、求解対象の問題を表すイジングモデル１に印加する磁場を低減させたときのイジングモデル１の状態変化のシミュレーションをする。処理部１２は、イジングモデル１の基底状態を求解する際に、シミュレーションで使用する係数の一部にノイズに相当する値を付加する。そして処理部１２は、イジングモデル１に印加する磁場の強さを、シミュレーション上の時間進行に伴って低減させていく実数時間発展処理（第１処理）を行う。実数時間発展処理によって、イジングモデル１に印加する磁場が徐々に低減され、その磁場が０になったとき、イジングモデル１が基底状態となる。

【0019】

また処理部１２は、実数時間発展処理の過程で、虚数時間発展法に基づき、イジングモデル１のエネルギーを低減させる熱的散逸処理（第２処理）を実行する。すなわち処理部１２は、シミュレーション上で進行させる時間を虚数化し、虚数時間を進行させることで、その時点で印加されている磁場に応じたイジングモデル１の基底状態を求める。磁場に応じたイジングモデル１の基底状態は、その時点の実数時間を固定パラメータとみなした定常状態におけるイジングモデルの基底状態（定常基底状態）である。実数時間発展処理が進行し、磁場の強さが０になったときのイジングモデル１の基底状態が、求解結果の基底状態である。

【0020】

このように初期状態にノイズを付加して実数時間発展と虚数時間発展とを組み合わせた計算を行うことにより、イジングモデル１の基底状態の求解を効率的に行うことができる。その理由について、以下に説明する。

【0021】

例えばＮ個（Ｎは自然数）の２値変数で定義された最小値求解問題において、量子力学に原理を置く計算を行い系の基底状態（最適値）および基底状態の波動関数（最適値における状態変数の組み合わせ）を求めるものとする。この方法において、ハミルトニアンとして、ＱＡ法におけるハミルトニアンを用いることができる。ＱＡ法では、初期状態はＮ個のスピンに横磁場をかけた状態である。例えば初期状態として実験的に準備可能な量子状態が用意される。そこから時間とともに横磁場を弱めていき、最終的に横磁場が完全になくなると、求解対象のハミルトニアンの基底状態、すなわち、エネルギー最小値が求まるようになっている。そして、数値的な解法においては標準的に量子モンテカルロ法（ＱＭＣ：Quantum Monte Carlo method）が使用される。ＱＡ計算では系の基底状態に保ち、励起状態に励起させないままゆっくり磁場を弱めることとなる。これは計算に対する非常に強い制約条件になっている。さらに１次の量子相転移と呼ばれる現象があり、磁場を弱める過程の中間状態で基底状態と励起状態のエネルギー差が指数関数的に小さくなる場合がある。この場合、横磁場の弱め方は指数関数的に遅くすることが要求される。

【0022】

この１次の量子相転移については近年、多体相互作用が１次の量子相転移を加速することが報告されている（特許文献１、非特許文献６参照）。
また、量子イジングハミルトニアンとスピン反転を記述する横磁場項で記述される時間依存シュレーディンガー方程式を虚数時間で解くことにより、量子系の動的挙動をシミュレーションできる。これは量子化されたシミュレーテッドアニーリングに対応するものであり、ＱＴＳＡ（Quantum Thermal Simulated Annealing）法と呼ぶこととする。ＱＴＳＡ法では、量子系の時間発展を直接シミュレーションするため、非平衡状態と平衡状態の区別はない。

【0023】

しかし、ＱＴＳＡ法をそのまま量子アニーリング計算に適用すると、基底状態を求めることができない場合がある。そこで、第１の実施の形態においては、最適化装置１０に以下の２つの概念を導入する。

【0024】

１つ目は量子試行である。量子アニーリングハミルトニアンにより波動関数の初期状態が時間発展をして終状態になると考えられ、一種の散乱問題と考えられる。最適化装置１０では、虚数時間における波動関数の散乱問題を考え、系の初期状態に依存して基底状態を含む励起状態を散乱問題として記述し、終状態を求める。

【0025】

２つ目は系の初期状態に多様さを持たせる技術である。これは通常ゲート型の量子コンピュータでも問題になるノイズを、組み合わせ最適化問題に積極的に利用することである。通常、ノイズは系の情報を破壊する有害なものであるが、ＱＴＳＡ法による計算で初期状態にノイズを乗せると、散乱問題の初期値に多様性を持たせることができる。その結果、散乱問題を解いたときに、有限確率で基底状態にたどり着くようにすることができる。

【0026】

ノイズの付加方法としては、例えば処理部１２は、イジングモデル１の初期状態に対してノイズに相当する値を付加し、ノイズに相当する値が付加されたイジングモデル１の状態をシミュレーション開始時の状態（ノイズ付き初期状態）とする。このとき処理部１２は、例えばイジングモデル１に含まれる複数のスピンのうちの一部のスピンの波動関数に含まれる係数にノイズに相当する値を付加する。なお処理部１２は、ノイズに相当する値を付加する一部のスピンを、例えば複数のスピンの中からランダムに選択する。

【0027】

処理部１２は、例えばイジングモデル１に対するノイズに相当する値の付加を複数回実行することで、異なるノイズが付加されたイジングモデル１の複数の状態を生成する。そして処理部１２は、イジングモデル１の複数の状態それぞれをシミュレーション開始時の状態として第１処理と第２処理とを実行する。すなわち、ノイズ付き初期状態の数だけ、シミュレーションが実行される。

【0028】

このようにイジングモデル１に様々なパターンのノイズを付加することで系の波動関数にノイズ付き初期状態に多様性を持たせることができる。ノイズ付き初期状態に多様性があることにより、一部のノイズ付き初期状態に対応する終状態が基底状態となる可能性が高くなる。その結果、ノイズが付加されていない初期状態を用いてシミュレーションを行った場合には終状態が基底状態とならないような場合であっても、イジングモデル１の基底状態を求めることが可能となる。またノイズの付加により系の対称性が破壊されていることで、量子アニーリングの計算における計算量の致命的な増大が緩和される。すなわち基底状態を容易に求めることができ、求解処理の効率化が図られている。

【0029】

また、厳密な意味での基底状態でなくとも、特定の少数個のスピンの０と１が異なっている状態が、準最適解として求められることも多い。そこで処理部１２は、さらに古典的なＭＣＭＣによる求解を行ってもよい。量子力学に基づく場合、スピンの状態は、スピンの取り得る状態（例えば０と１）の量子力学的重ね合わせによって表現される。量子力学的重ね合わせで表現されたスピンの状態は、観測した場合に各状態が観測される確率を示している。それに対して、古典力学に基づくイジングモデル（古典イジングモデル）では、スピンの状態は、取り得るいずれかの状態となる。

【0030】

例えば処理部１２は、スピンの取り得る状態の量子力学的重ね合わせで各スピンの状態が表現されたイジングモデル１に対して、量子力学に基づいて第１処理と第２処理とを繰り返し実行する。次に処理部１２は、イジングモデル１のエネルギーが収束した後のイジングモデル１の状態に基づいて、スピンの状態を取り得る値のいずれかに決定した古典イジングモデルの状態を求める。そして処理部１２は、古典イジングモデルの状態を用いた古典力学に基づく解法（例えばＭＣＭＣ）によって古典イジングモデルの基底状態を求める処理（第３処理）を実行する。

【0031】

このように古典力学に基づく解法と組み合わせることで容易に基底状態にたどり着くことができる。すなわち、第１処理と第２処理だけでは厳密な基底状態が求められない場合でも、例えば古典力学に基づくＭＣＭＣ（以下、古典ＭＣＭＣと呼ぶこともある）では高々数ステップで基底状態にたどり着く状態の場合がある。このような場合、古典ＭＣＭＣと組み合わせることで、基底状態を求めることが可能となり、求解効率を向上させることができる。

【0032】

系の対称性を破壊して計算を加速するもう一つの方法として横磁場項の強さにノイズを入れる方法もある。横磁場の強さにノイズをいれることでも波動関数の対称性を破壊できる。例えば処理部１２は、イジングモデル１に含まれるスピンにおける磁場の強さを示す係数にノイズに相当する値を付加する。これにより系の対称性が破壊され、量子アニーリングの計算における計算量の致命的な増大が緩和される。

【0033】

〔第２の実施の形態〕
第２の実施の形態は、ノイマン型コンピュータを用いて、２値Ｎ変数で定義された最小値求解問題について、量子力学に原理を置く計算を行い系の基底状態（最適値）および基底状態の波動関数（最適値における状態変数の組み合わせ）を求めるものである。

【0034】

＜最適化装置のハードウェア構成＞
図２は、最適化装置のハードウェアの一例を示す図である。最適化装置１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。

【0035】

メモリ１０２は、最適化装置１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に利用する各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。

【0036】

バス１０９に接続されている周辺機器としては、ストレージ装置１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。

【0037】

ストレージ装置１０３は、内蔵した記録媒体に対して、電気的または磁気的にデータの書き込みおよび読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。

【0038】

グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、有機ＥＬ（Electro Luminescence）を用いた表示装置や液晶表示装置などがある。

【0039】

入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。

【0040】

光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取り、または光ディスク２４へのデータの書き込みを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ－ＲＡＭ、ＣＤ－ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ－Ｒ（Recordable）／ＲＷ（ReWritable）などがある。

【0041】

機器接続インタフェース１０７は、最適化装置１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。

【0042】

ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。

【0043】

最適化装置１００は、以上のようなハードウェアによって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示した最適化装置１０も、図２に示した最適化装置１００と同様のハードウェアにより実現することができる。

【0044】

最適化装置１００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。最適化装置１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、最適化装置１００に実行させるプログラムをストレージ装置１０３に格納しておくことができる。プロセッサ１０１は、ストレージ装置１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。また最適化装置１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばプロセッサ１０１からの制御により、ストレージ装置１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。

【0045】

このようなハードウェアの最適化装置１００を用いて、２値Ｎ変数で定義された最小値求解問題について、量子力学に原理を置く計算を行い、系の基底状態および励起状態の波動関数を求めることができる。ただし、前述のように、ＱＡ法を用いた解探索の時間が長期化する原因の１つとして、１次の量子相転移と呼ばれる現象の発生がある。すなわち１次の量子相転移の発生を抑制しようとすると、解探索の時間が長期化してしまう。そこで、１次の量子相転移の発生を抑止しようとすると、最小値求解問題の解探索時間が長期化する理由について、詳細に説明する。

【0046】

なお、イジングモデルの基底状態とは、イジングモデルを表すハミルトニアンのエネルギーが最も低い状態である。以下、イジングモデルの基底状態を、ハミルトニアンの基底状態と呼ぶこともある。また、ＱＡ法においてイジングモデルに印加される磁場は横磁場と呼ばれる。

【0047】

＜基底状態の算出手法の概要＞
ＱＡ法で求解する場合、横磁場を加えた項が初期状態として用意される。熱的散逸処理を行わない場合、ここから、求めたい２値Ｎ変数問題として定義された最適化問題のハミルトニアンの横磁場の強さを十分ゆっくり弱めていき、最終的に横磁場が０になるようにする。このとき、系の波動関数を初期状態から求解対象のイジングモデルの基底状態の波動関数に十分ゆっくり遷移させていき、励起状態に励起させないように磁場を弱めることとなる。これを断熱遷移と言う。断熱遷移の結果として、最終的に系は求解対象のハミルトニアンの基底状態に遷移する。このようにして求解対象の量子イジングハミルトニアンの基底状態を求めることができる。

【0048】

一方でＱＡ法では、系のハミルトニアンによっては基底状態と励起状態のエネルギーギャップが小さくなり、ハミルトニアンが励起状態に遷移しないようにするため、断熱遷移にかける時間が伸びてしまうという問題がある。すなわち、励起状態への遷移は、断熱遷移をしていく過程で、中間のハミルトニアンにおける基底状態と第一励起状態のエネルギー差が小さくなってしまうことに起因して発生する。ＱＡ法を用いた場合、励起状態への遷移を抑止するために、最適化装置１０は、基底状態と第一励起状態のエネルギー差が小さいときには、励起状態への遷移が発生しないように、横磁場を弱める量を少なくして、ハミルトニアンをゆっくり遷移させることとなる。

【0049】

ここで、エネルギー差がスピンの数の関数として指数関数的に小さくなる場合を１次の量子相転移と呼ぶ。量子相転移が１次の場合、ＱＡ法による断熱遷移条件を守りながら時間発展するための時間刻み幅は、基底状態と励起状態のエネルギーギャップに依存する。そのため最適化装置は、時間刻み幅を指数関数的に小さくすることとなる。これはＱＡ法により終状態に到達するための時間発展回数が指数関数的に延びることを意味する。

【0050】

他方において、エネルギー差がスピンの数の関数としてべき乗で小さくなる場合も存在する。これを２次の量子相転移と呼ぶ。２次の量子相転移が起こる場合、ＱＡ法の断熱遷移の速度はべき乗で小さくなる。この場合は終状態にたどり着くために要する時間は多項式時間で増加する。

【0051】

しかし、応用上興味のある問題の多くが１次の量子相転移を伴うため、多くの実用上の問題において計算時間が指数関数的に増大する。そこで１次の量子相転移を２次の量子相転移に置き換えることが可能であるかを探る研究もなされている。Non-Stoquastic Hamiltonianと呼ばれる特殊なハミルトニアンの系において、横磁場項に工夫をすることで１次相転移を２次相転移にすることができたという報告もある（非特許文献６）。従って、通常のＱＡ法の枠組みにおいては１次の量子相転移をいかにして回避するかが計算を実行する上での中心的な課題になる。

【0052】

さらにＱＡにおいては断熱遷移条件を守ることが求められている。しかし、ＱＡ法に基づく量子コンピュータを用いた計算では、この断熱条件が破られているにも拘らず、有限の確率で基底状態を計算できている可能性があることが報告されている（非特許文献７）。量子アニーリングの実機を用いた計算においては、最適なＴＴＳ（Time To Solution）が存在するという報告であるが、物理学的な見地において、なぜ最適なＴＴＳが存在するかについては一考の余地がある。

【0053】

これらの状況を考えると、量子アニーリングにおいては、何らかの原因でエネルギー散逸（熱的散逸）が起こり、エネルギー散逸によって量子アニーリングの機構が機能している可能性があると考えられる。従って、量子論の枠組みにおいて熱的散逸が起こりうるものとして取り扱うことが、ＱＡプロセスを本質的に理解するために重要となる。

【0054】

熱的散逸は外部環境との相互作用が本質である。そのため、ハミルトニアンに、外部環境の影響を熱的散逸項で取り込むことが求められる。原子物理学分野ではこのエネルギー散逸を取り扱う方法が開発されている。それは時間依存シュレーディンガー方程式を虚数時間で解くことである。虚数時間を導入することで波動方程式が拡散方程式になることが、エネルギー散逸を取り扱うことが可能となる主な理由である。

【0055】

エネルギー散逸を取り扱う方法の具体例を挙げると、稀薄原子気体のボース・アインシュタイン凝縮（ＢＥＣ：Bose Einstein Condensation）によるボース凝縮体の取り扱いのために考えられた方法がある。稀薄原子気体のＢＥＣの理論は多面的に研究されているが、理論的な取り扱いで厄介になるのは量子渦の計算である。そして、ボース凝縮体を実験的に光学的なスプーンでかき混ぜると、量子化された渦が実験的に観測されることが分かっている（非特許文献８）。この量子渦のボース凝縮体の基底状態は平均場近似の結果として得られるＧＰ（Gross-Pitaevskii）方程式に角運動量項を入れた方程式でよく近似されることが分かっている。

【0056】

しかし、渦格子のＧＰ方程式の基底状態を求めることは非常に難しい。なぜならば、渦格子の対称性も渦の個数も方程式からは何もわからないからである。このように基底状態が存在するにも関わらず、対称性などの情報が一切わからない状況が存在する。このときに渦格子の基底状態を計算するために使用される方法が虚数時間発展法である。虚数時間発展法を用いることで、系の基底状態に対して対称性などの情報がない場合でも基底状態を計算できる（非特許文献９）。また、この方法は現象論的な熱散逸項を導入するためにも使うことができ、ＢＥＣ中のDark-Soliton状態を扱った研究（非特許文献１０）などもある。虚数時間発展法は簡単であり、かつ、強力なアルゴリズムであるが、量子アニーリングの系では運動エネルギー演算子が明示的に表れないため、量子アニーリングの系に虚数時間発展法を適用するのは容易ではない。

【0057】

虚数時間発展法においては、時間を虚数にすることで、量子力学的な時間発展演算子がｅｘｐ（－ｉＥｔ）型からｅｘｐ（－Ｅτ）型に変わる（ｉは虚数単位、Ｅはエネルギー、ｔは実数時間、τは虚数時間である）。このため、任意の波動関数は完全系｜Ψ_n＞（ｎ番目の固有状態）を用いて式（１）のように展開できるが、時間を虚数にすることで式（２）のようになる。

【0058】

【数1】

【0059】

【数2】

【0060】

｜Ψ（ｔ）〉は時間ｔのスピンの状態を表す波動関数であり、ｄ_n（ｔ）は｜Ψ_n〉の係数であり、ｉは虚数単位であり、Ｅ_nはハミルトニアンのｎ番目のエネルギーであり、Ｅ₀は基底状態のエネルギーであり、ｈバーはディラック定数であり、τは虚数時間である。式（２）において、虚数時間τ→∞の極限を取ると、基底状態のみの振幅が残り、式（３）のようになる。

【0061】

【数3】

【0062】

ただし、このままでは基底状態の振幅もｅｘｐ（－Ｅ₀τ／ｈ）（ｈはバー付き）に従い０になってしまう。これは時間を虚数にすることでユニタリ性（確率の保存）が破れるためである。そのため、計算の都度波動関数の規格化をやり直し、ユニタリ性を強制的に守るようにして計算を継続する。このようにすると基底状態にたどり着くことができる。

【0063】

この方法は変分原理に基礎をおいており、問題に依存せず適用することができる。変分原理に基礎をおいているため、最適解に到達することが原理的には保証されている方法である。虚数時間発展法で取り扱うことにより、量子イジング系の時間発展を決定論的に取り扱うことができる。

【0064】

図３は、虚数時間発展法を用いた熱的散逸の模式図である。図３には、横軸にλ、縦軸にハミルトニアンによって求められる系のエネルギーをとったグラフ３０が示されている。図３の例では、実数時間発展に伴い横磁場の強さに対応する定数λが０から１に遷移するときの、λに応じた基底状態におけるエネルギーを実線３１で示している。また、λに応じた第１励起エネルギー状態を点線３２で示している。λが０の時に横磁場の強さが最大であり、λが１の時に横磁場の強さが最小である。

【0065】

通常のＱＡ法では磁場を大きく変えてしまうと、点Ｐで基底状態にあった状態が、第１励起エネルギー状態の点Ｑに遷移してしまう。通常のＱＡ法では、このような状況を許さないため、磁場を断熱的に遅く緩めることとなる。しかし、点Ｑからエネルギーを脱励起させて点Ｒに落とすことができるのであれば、断熱遷移条件は厳密に守ることは要求されない。最適化装置１００は、虚数時間発展により、点Ｑからエネルギーを脱励起させて点Ｒに落とすことを可能とする。

【0066】

一方で、虚数時間発展法を用いたとしても、１次の量子相転移の問題は起こる。１次の量子相転移が起こると、量子アニーリング中の横磁場を弱めた場合の中間状態において、基底状態と励起状態のエネルギー差がほぼ縮退した状態になってしまう状況が生じる。基底状態と励起状態のエネルギーが縮退すると、エネルギー的には状態が区別できなくなる。結果として、基底状態に到達するか、励起状態に到達するかを制御不能になってしまう。すなわち、虚数時間発展を用いたとしても、必ずしも基底状態が求められるわけではなく、一定確率で解が求まることになる。これは系の波動関数の初期状態の問題でもあるが、系のハミルトニアン自体にも問題がある。

【0067】

そのため、量子アニーリング法において、ＱＴＳＡの枠組みで基底状態に到達するための技術開発には、基礎理学としての価値と工学としての価値の両方がある。第２の実施の形態は、基底状態を系統的に求めるための技術を発展させたものである。

【0068】

＜イジングモデルと問題のハミルトニアン＞
次に、最適化装置１００による求解の対象となる、イジング型の最小値求解問題（イジング型問題）について説明する。イジング型問題は元々物理学領域において磁性体の研究に使われたモデルであり、Ｎ個のスピンがあるときの全エネルギーが次の式（１．１）で記述される。

【0069】

【数4】

【0070】

右辺の１項目は、Ｎ個の状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘ_iはｉ番目の状態変数、ｘ_jはｊ番目の状態変数を表し、Ｗijは、ｘ_i，ｘ_jの相互作用の大きさを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂ_i）と各状態変数の値との積の総和を求めたものであり、右辺の３項目（Ｃ）は定数である。

【0071】

最小値求解問題は、式（１．１）で与えられるハミルトニアンの最小値を求める問題である。なお、式（１．１）に示すハミルトニアンの最小値をモンテカルロ法などに基礎を置く計算で求めると、最小値を求めるための計算量が、一般に２^N回程度になることが知られている。つまり、最適解を求めるための計算時間は多項式時間で抑えられず、指数関数的に増大する。

【0072】

この問題を解決するため量子論が導入される。式（１．１）に示すハミルトニアンを量子化されたハミルトニアンに直したものが次の式（１．２）、式（１．３）、式（１．４）である。

【0073】

【数5】

【0074】

【数6】

【0075】

【数7】

【0076】

式（１．２）は式（１．１）を量子化したものである。式（１．２）の左辺のハミルトニアンＨ_TH（添え字のＴＨはTarget Hamiltonianの略）である。式（１．３）、式（１．４）のσ_i ^z（σは＾付き）はｉ番目のスピンのパウリ行列のｚ成分である。

【0077】

なお、式（１．２）、式（１．３）、式（１．４）だけでは量子力学のダイナミクスを記述したときにスピン反転を記述することができない。そこで横磁場項を加えた式（１．５）のハミルトニアンを量子力学的なハミルトニアンとして定義する。

【0078】

【数8】

【0079】

式（１．５）におけるＴ（＾付き）は横磁場項と呼ばれ、時間依存した量子力学の方程式に直したときにスピン反転を記述するために導入された項であり、次の式（１．６）で定義される。

【0080】

【数9】

【0081】

Ｊは、スピンの強度を規定する係数である。またハミルトニアン中のＧ₁は縦磁場項の強さを規定する定数、Ｇ₂は横磁場項の強さを規定する定数である。式（１．６）のσ_i ^x（σは＾付き）に関して１次の項のみを取り出した場合が非特許文献５によって導出されたＱＡ法である。非特許文献５のＱＡ法ではＧ₁＝λ、Ｇ₂＝１－λにとる。λは、０≦λ≦１の横磁場の強さを示すパラメータである。一般にＧ₁とＧ₂の間に関係性を持たせずに、値を独立にとってよい。

【0082】

次に問題の演算子の行列、ベクトル表記を以下の式で定義する。

【0083】

【数10】

【0084】

【数11】

【0085】

【数12】

【0086】

【数13】

【0087】

式（１．７）は通常用いられる定義と異なっているが、これはスピン演算子の固有値を１および－１ではなく、１および０にするためである。｜α_i＞と｜β_i＞には直交関係として以下の式が成立する。

【0088】

【数14】

【0089】

次に縦磁場項σ_i ^z（σは＾付き）についての固有方程式を以下に示す。

【0090】

【数15】

【0091】

【数16】

【0092】

横磁場項σ_i ^x（σは＾付き）については以下のようになる。

【0093】

【数17】

【0094】

【数18】

【0095】

次にｉ番目のＱｂｉｔを表す１粒子波動関数のケット表示を式（１．１６）で定義する。

【0096】

【数19】

【0097】

Ｃ_iα（ｔ）は、ｉ番目のスピンの｜α_i〉の状態の確率振幅と位相を表す複素数である。Ｃ_iβ（ｔ）は、ｉ番目のスピンの｜β_i〉の状態の確率振幅と位相を表す複素数である。１粒子波動関数のブラ表示を次の式（１．１７）で定義する。

【0098】

【数20】

【0099】

Ｃ_iα ^*（ｔ）、Ｃ_iβ ^*（ｔ）はそれぞれＣ_iα（ｔ）、Ｃ_iβ（ｔ）の複素共役である。量子力学において波動関数は一般に複素数で定義された関数であるため、係数も複素数として定義される。ここで、１粒子波動関数は規格化されているとして、式（１．１８）が成り立つ。

【0100】

【数21】

【0101】

次によく使用する物理量を以下のように定義する。

【0102】

【数22】

【0103】

【数23】

【0104】

Ｎ粒子系の波動関数を以下のように定義する。

【0105】

【数24】

【0106】

なお、ｍ番目のスピンがαスピンとして観測される確率をＰ_m,α、βスピンとして観測される確率をＰ_m,βとする。Ｐ_m,αとＰ_m,βとは以下の通りである。

【0107】

【数25】

【0108】

【数26】

【0109】

＜虚数時間発展法による基礎方程式＞
今、ハミルトニアンの式（１．５）で定義される時間依存シュレーディンガー方程式である式（２．１）を解くことを考える。

【0110】

【数27】

【0111】

式（２．１）のｎ番目の固有値をＥ_n、固有状態｜Ψ_n＞は次の式（２．２）の解である。

【0112】

【数28】

【0113】

任意の波動関数｜Ψ（ｔ）＞は完全系｜Ψ_n＞で以下のように展開できる。

【0114】

【数29】

【0115】

ここで、時間を以下の式（２．４）を用いて実数から虚数時間τに変換する。

【0116】

【数30】

【0117】

虚数時間に対して、式（２．３）は次の式（２．５）になる。

【0118】

【数31】

【0119】

ここで、虚数時間τ→∞の極限を取ると、励起状態の成分は基底状態よりも高速に減衰するため、式（２．６）のように基底状態の成分のみが残ることが分かる。

【0120】

【数32】

【0121】

しかし、係数因子ｅｘｐ（－Ｅ₀τ）ｈ）（ｈはバー付き）自体もτ→∞の極限で０になる。これは時間を虚数にすることで、時間発展演算子がユニタリではなくなるためである。具体的には波動関数の規格化積分が１以下になってしまうためである。このため、計算の進行と共に強制的に規格化積分を１にするようにしてしまう。この操作を再規格化（Renormalization）と呼ぶ。虚数時間発展と再規格化を用いることで系の基底状態にたどり着くことが原理的には保証される。

【0122】

次に虚数時間発展法（ＩＴＰ：Imaginary Time Propagation）を量子イジング系に適用した場合の支配方程式を導出する。時間発展演算子として１次のＩＴＰを用いる。１次のＩＴＰの演算子は式（２．７）で与えられる。

【0123】

【数33】

【0124】

式（２．７）では、Δτの２次の項は関数Ｏ（Δτ²）で纏めている。式（２．７）をΔｔに関して１次の項まで展開して整理すると式（２．８）となる。

【0125】

【数34】

【0126】

次に波動関数の時間発展を記述する方程式を導く。波動関数の時間発展は式（２．９）で与えられる。

【0127】

【数35】

【0128】

今、ｍ番目の粒子の波動関数だけ除いたＮ－１粒子系の波動関数を｜Ψ_m（τ）＞（ｍはオーバーライン付き）と表記すると、次の式（２．１０）のようになる。

【0129】

【数36】

【0130】

式（２．９）の両辺に＜α_m｜×＜Ψ_m（τ₀）（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き）を作用させる。その結果、式（２．９）の左辺は式（２．１１）となる。

【0131】

【数37】

【0132】

同様に式（２．１０）の両辺に＜β_m｜×＜Ψ_m（τ₀）（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き）を作用させる。その結果、式（２．９）の左辺は式（２．１２）となる。

【0133】

【数38】

【0134】

ここで、Ｍ_i ^*はスピンの変化のしにくさを表す量であり、以下の式（２．１３）で表される。

【0135】

【数39】

【0136】

ここで以下の式（２．１４）を定義する。

【0137】

【数40】

【0138】

Ω_m（τ）をΔτの項まで整理すると式（２．１５）となる。

【0139】

【数41】

【0140】

次に式（２．９）の両辺に＜α_m｜×＜Ψ_m（τ₀）と＜β_m｜×＜Ψ_m（τ₀）（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き）を作用させた結果の右辺を求め、支配方程式を求めると、式（２．１６）となる。

【0141】

【数42】

【0142】

ここで、行列要素は以下の式（２．１７）～（２．１９）で与えられる。

【0143】

【数43】

【0144】

【数44】

【0145】

【数45】

【0146】

式（２．１７）～式（２．１９）中のパラメータＥ_TH，ｈ_z,m，Ｔ_αα，Ｔ_ββ，Ｔ_αβを以下に与える。
Ｅ_THは求解対象のハミルトニアンのエネルギー項であり、式（２．２０）～式（２．２２）で定義される。

【0147】

【数46】

【0148】

【数47】

【0149】

【数48】

【0150】

またｈ_z,mは平均場に相当する量であり、式（２．２３）、式（２．２４）で与えられる。

【0151】

【数49】

【0152】

【数50】

【0153】

Ｔ_ααは横磁場項の期待値であり、＜α_m｜×＜Ψ_m（τ₀）｜Ｔ｜Ψ＞（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き、Ｔは＾付き）からＣ_mα（τ）の係数を抜き出したものである。Ｔ_αβは＜α_m｜×＜Ψ_m（τ₀）｜Ｔ｜Ψ＞（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き、Ｔは＾付き）からＣ_mβ（τ）の係数を抜き出したもの。そして、Ｔ_ββは＜β_m｜×＜Ψ_m（τ₀）｜Ｔ｜Ψ＞（×は丸の中に×の記号（テンソル積）、Ψの添字ｍはオーバーライン付き、Ｔは＾付き）からＣ_mβ（τ）の係数を抜き出したものである。

【0154】

行列形式の評価は次のようにして行う。ここで、ｍ番目の１粒子波動関数だけが｜φ_m＞に変わったＮ粒子波動関数として式（２．２５）を導入する。

【0155】

【数51】

【0156】

ここで、｜φ_m＞を式（２．２６）で定義する。

【0157】

【数52】

【0158】

いま、式（２．２５）から次の積分を考える。

【0159】

【数53】

【0160】

式（２．２６）の積分は式（２．２７）のように評価することができる。

【0161】

【数54】

【0162】

このようにして評価した行列形式が、式（２．１６）における波動関数の時間発展の右辺における２行２列の行列要素を与える。いまＴ_αα，Ｔ_ββ，Ｔ_αβの行列要素は、式（２．２９）、式（２．３０）として与えられる。

【0163】

【数55】

【0164】

【数56】

【0165】

なお行列要素Ｔ_αα，Ｔ_αβ，Ｔ_βα，Ｔ_ββは積分＜φ（ｔ）｜Ｔ_n｜Ψ（ｔ）＞（Ｔは＾付き）から計算できる。
いま行列要素を具体的に与えると式（２．３１）となる。

【0166】

【数57】

【0167】

ここで、Ｘ_mは式（２．３２）とする。

【0168】

【数58】

【0169】

次にΩ_m（τ）について説明する。αスピンが観測される確率（１．２２）、βスピンが観測される確率は（１．２３）で与えられる。従って、観測確率の和を虚数時間τについて微分すると式（２．３３）となる。

【0170】

【数59】

【0171】

観測確率の和については、常にＰ_m,α＋Ｐ_m,β＝１であることを用いている。従って、Ｍ_m ^*は実数部を持たない純虚数である。
ここで、時間発展の支配方程式（２．１６）において、波動関数の初期値としてＣ_mα（τ＝０）とＣ_mβ（τ＝０）を共に実数に選んだ場合を考える。式（２．１６）の右辺の行列要素はすべて実数である。従って、Ｃ_mαとＣ_mβが実数であればＣ_mα（τ＋Δτ）Ω_m（τ）は実数である。一方で、Ω_m（τ）はＭ_m ^*から計算されるが、Ｍ_m ^*の定義式である式（２．２３）から時刻τにおけるＭ_m ^*も実数になる。ところが、式（２．１６）からＭ_m ^*は純虚数になるため、この場合はＭ_m ^*＝０になる。従って、Ω_m（τ）＝１になり、Ｃ_mα（τ＋Δτ）も実数になる。Ｃ_mβ（τ＋Δτ）もまったく同様の議論で実数になる。

【0172】

次に波動関数の初期値Ｃ_mα（τ＝０）、Ｃ_mβ（τ＝０）を複素数として決めた場合を考察する。この場合、Ｍ_m ^*は複素数として与えられる。この場合、一般には、Ω_m (τ)≠１であり、複素数になる。そして、Ｍ_m ^*はｍ番目のスピンから決まる量である。従って、Ω_m (τ)は着目するｍ番目以外の全てのスピンの情報から決まる。従って、Ｃ_mα（τ＋Δτ）は直前のｉ≠ｍを満たす全ての他のスピンの波動関数Ｃ_iα（τ）とＣ_iβ（τ）の影響を受ける。すなわち、時刻τ＝０の波動関数の情報に位相情報を組み入れ、複素数にした場合、虚数部分を介して情報をやり取りする。絶対値の２乗で定義される観測確率には影響を与えない位相部分がスピン間の相互作用にも重要な役割を果たしている。そして、この相互作用は他の全てのスピンと結合し、位相部分を介しても結合する。

【0173】

このようにして量子アニーリングを決定論的に導く支配方程式を導出できる。なお、全エネルギーについては式（２．３４）、式（２．３５）で与えられる。

【0174】

【数60】

【0175】

【数61】

【0176】

計算量については虚数時間発展の支配方程式（２．１６）を解くための計算量はＯ（２Ｎ）のオーダーとなる。式（２．１６）中の行列要素の計算には求解ハミルトニアンの全エネルギーの計算が必要であるため、全エネルギーの計算のための計算量はＯ（Ｎ²）のオーダーとなる。従って、虚数時間発展法の収束までの計算回数をＬとすると、虚数時間発展法ではＯ（ＬＮ²）の計算量が要求される。しかし、虚数時間発展においてエネルギー減衰はｅｘｐ（ΔＥ₀₁τ）で減衰収束すると評価できるため、指数関数的に高速にエネルギーは下がる。

【0177】

＜量子試行の定義＞
以下、量子試行という用語を用いる。量子試行は以下のように定義する。時刻ｔ＝０におけるＮ体波動関数｜Ψ₀＞に対して、虚数時間における時間発展演算子を式（３．１）とする。

【0178】

【数62】

【0179】

このとき、虚時間発展の結果を得られる終状態を｜Ψ（τ）＞とし、式（３．２）の関係が得られる。

【0180】

【数63】

【0181】

このとき、｜Ψ₀＞が異なれば、終状態も一般には異なる。初期状態から終状態を得る課程を量子試行と呼ぶことにする。このようにして、量子アニーリングを一種の散乱過程として捉えることが可能である。

【0182】

＜波動関数に対するｋ粒子励起的ノイズ付加による初期状態の生成＞
初期状態を複数の方法で用意できる。Ｎ粒子波動関数を構成する１粒子波動関数の初期値を以下のように準備する。

【0183】

【数64】

【0184】

ここで、式（４．１）ですべての１粒子波動関数の係数を同じにする初期値を使うものとする。すなわち、ｉ≠ｊなる任意の２スピンに対して、Ｃ_iα（ｔ＝０）＝Ｃ_jα（ｔ＝０）、Ｃ_iβ（ｔ＝０）＝Ｃ_jβ（ｔ＝０）になるような初期状態を準備する。このとき、ｋ粒子励起的ノイズを次のように定義する。Ｎ個のスピンの中からｋ個のスピンを選択し、そのそれぞれに対して、ｚパーセントのノイズを確率振幅に対して入れる。すなわち、１番目のスピンが選択されたとき、ノイズ付加後の確率振幅を式（４．２）および式（４．３）で定義する。

【0185】

【数65】

【0186】

【数66】

【0187】

式（４．２）においてＣ_iα＝０なる初期状態を用いた場合、ノイズが加わらないため、この場合式（４．２）、式（４．３）はＣ_iβに対してノイズを加えるものとする。
ただし、式（４．２）と式（４．３）によるノイズ付加の定義式は一例である。ノイズを付加することにより、他の量子スピンと特定のスピンだけが異なる１粒子量子状態を取ることが本質的に重要である。また、ｚの値は一般にはスピン毎に異なるようにしておく。なぜならば、ｋ個のスピンにノイズが加わるとしても、ｋ個すべてに同じノイズが加わる必然性が無いためである。そのため最適化装置１００は、乱数などを用いてノイズの振幅であるｚを決める。

【0188】

式（４．２）、式（４．３）式によってｋ個のスピンに対してノイズを付加した場合、確率振幅が増大または減少するため、他の状態が部分的に入ってきていると考えれば、これは励起を部分的に取り入れたと考えることができる。そのため、これをｋ粒子励起的ノイズと呼ぶことにする。なお、ノイズ無しの初期状態としてＣ_iα（ｔ＝０）≠Ｃ_jα（ｔ＝０）、Ｃ_iβ（ｔ＝０）≠Ｃ_jβ（ｔ＝０）になるように、最初から対称性を崩した初期状態を準備しても良い。例えば、一つのスピンについてＣ_iα（ｔ＝０）=１にし、残りのスピンに関してＣ_jα（ｔ＝０）=０（ただし、ｉ≠ｊ）として選んでもよい。このように選んだ場合でも、ノイズを付加すると終状態は異なるからである。

【0189】

＜量子試行とノイズ付加により解決を目指す課題について＞
ここまでの議論においては、量子試行と波動関数に対するノイズを天下り的に与えた。この２つを導入することで、虚数時間発展だけでは解決されていない課題が解決されるメカニズムを説明する。虚数時間発展だけでは解決されていない課題とは、虚数時間発展の枠組みにおいて量子系の時間発展をしても、必ずしも基底状態を求めることができないことである。

【0190】

ここで、虚数時間発展で実際に４都市の巡回セールスマン問題（ＴＳＰ：Traveling Salesman Problem）を解いたときに生じる具体的な問題点を提示する。
図４は、４都市のＴＳＰにおける都市の配置例を示す図である。図４の例では、１番目の都市の位置が（０，０）、２番目の都市の位置が（３，０）、３番目の都市の位置が（３，１）、４番目の都市の位置が（１，１）である。

【0191】

図４に示した問題の場合、解の候補は３通りある。各候補の経路長は、短い順から7．４１４、７．８１２、および１０．３９８である。
ＴＳＰをイジングモデルに変換し、非特許文献５にあるような外部磁場の下で量子アニーリング計算を、虚数時間発展法を用いて行うこととする。なおＧ₁＝１－Ｇ₂とし、０≦Ｇ₂≦１とする。虚数時間τ＝０でＧ₂＝１とする。Ｇ₂は離散的な点として与え、等間隔に分点を定義する。ここで、計算に用いる分点の数ｎを６４から１０２４まで増やして計算を行う。ｎ番目のＧ₂値をＧ₂（ｎ）と書く。中間状態は式（５．１）で定義する。

【0192】

【数67】

【0193】

式（５．１）の中間状態のそれぞれについて、虚数時間発展を実行し、エネルギーが収束するまで計算する。計算収束の閾値は１０^-10とする。波動関数の初期値はこの計算では任意であるが、Ｃ_mα（τ＝０）＝１、Ｃ_mα（τ＝０）＝０として初期化する。また、Ｇ₂（０）＝１、Ｇ₁=１．０とした。ｎ＝０での計算が終了すると、ｎの値を１増やし、Ｇ₂の値が減少した次の中間状態Ｇ₂（１）、Ｇ₁（１）での基底状態を求めるために虚数時間発展を行う。このとき、新たな中間状態での波動関数の初期状態は前段の最終結果として得られたＧ₂（０）においてエネルギーが収束した波動関数を用いる。この計算を繰り返し、最終的にＧ₂（ｎ＝１００）＝０の計算が終了すると、求めるべき求解ハミルトニアンの基底状態（または励起状態）が求まる。この一連の計算が量子試行である。

【0194】

図５は、Ｇ₂の値の数に応じた計算結果の違いを示す図である。図５に示すグラフ３３は横軸がＧ₂であり、縦軸がエネルギーである。グラフ３３に示す３本の実線それぞれが、異なる分点数の量子試行によるＧ₂の値ごとのエネルギー値を示している。グラフ３３に示すように、分点の数によっては結果が７．４１４の基底状態が求まる場合もあれば、７．８１２や１０．３９８の励起状態が求まる場合もある。

【0195】

４都市ＴＳＰの場合は打切り誤差を小さくし、分点数を多くすれば確実に基底状態を求めることは可能である。しかし、ＴＳＰＬＩＢ（ＴＳＰの問題集）に収録されている問題のうちの実用上の問題に近いulysses16などを解くと状況は悪化する。エネルギーの打切り誤差を倍精度浮動小数点の桁数にとり、分点数をかなり大きく取っても基底状態に行き着かないのである。

【0196】

このように基底状態に行き着かない可能性があるものの、４都市のＴＳＰで求める３つの解候補はすべて計算できている。この計算は断熱遷移条件を守っていると言い難いが、一定確率で基底状態が求まるという量子アニーリングの実験結果を支持する。そのため、量子アニーリング計算は実機では熱的散逸による自然なエネルギー散逸から計算が行われていると解釈することができる。

【0197】

このような状況が生じるメカニズムを次に説明する。
図６は、横磁場の強さに応じたエネルギー値を模式的に示す図である。図６に示すグラフ３４は、横軸が横磁場の強さであり、縦軸はエネルギーである。グラフ３４には典型的な１次の量子相転移と呼ばれる状況が表されている。すなわち、量子アニーリング時に横磁場を弱めて行くと、求解対象のハミルトニアンの影響が強くなっていき、複数の量子状態が実質的に縮退してしまう現象が表されている。定性的にはエネルギー差がΔＥ≒ｅｘｐ（－ｂＮ）のように指数関数的に小さくなってしまう（Ｎはスピン数、ｂは係数である）。

【0198】

量子アニーリングはグラフ３４中の一番下のポテンシャル曲線で表される量子状態を保ちながら横磁場を弱めることで求解対象の基底状態にたどり着く計算である。しかし、量子ビット数（スピン数）が多くなると、エネルギー差は指数関数的に小さくなる。そのため、エネルギー散逸が起こり外界と相互作用をしたとしても、エネルギー差が実質上０になる状況が存在してしまう。数値計算の場合、倍精度浮動小数点がよく使用されるが、これは１６桁しか桁数がない。そのため、１７桁目以降に違いがでるような状況が存在すれば、倍精度浮動小数点ではそのような状況は記述できない。すなわち、数値的に全て縮退した状態として扱われてしまう。

【0199】

そこで第２の実施の形態では発想を逆転させ、エネルギーが分解能の範囲で縮退しているならば、縮退したどの状態も等しく実現されると考える。つまり、１００％の確率で求めるのはあきらめ、有限確率で求めることができればよいと考えるのである。そして、このように考えれば、励起状態も求められると考えるのである。しかし、計算手続きが同じであれば同じ微視的状態に落ち着くと考えられるため、意図的に微小なノイズを付加して初期状態に変調を与える。ＱＴＳＡの支配方程式である式（２．１６）はポテンシャル項に現時刻での自身の波動関数を含む一種の平均場近似の方程式になっているため、典型的な複雑系としてふるまう。そのため、ノイズをわずかに加えただけでも、そのノイズは時間発展中に増幅され系の予期しない振る舞いを引き起こす。系は確率的に求解対象の基底状態にも励起状態にも行き着くのである。第２の実施の形態に係る最適化装置１００は、系は確率的に求解対象の基底状態にも励起状態にも行き着くようにするために量子試行とノイズを導入する。

【0200】

図７は、ノイズ付加を伴うＱＴＳＡ計算を模式的に示す図である。例えば最適化装置１００は、まずスピンの状態を示すスピン情報４０として所定の初期値を設定する。スピン情報４０では、例えば各スピンの向きが０または１に設定されている。最適化装置１００はスピン情報４０に対してノイズを付加し、複数のスピン情報４１，４２，４３，・・・を生成する。スピン情報４１，４２，４３，・・・それぞれは、一部のスピンがランダムに選択され、選択されたスピンについて０と１以外の状態（量子力学に基づく確率表現）に変更されている。

【0201】

最適化装置１００は、スピン情報４１，４２，４３，・・・それぞれを初期状態として、ＱＴＳＡの計算を行う。なおスピン情報４１，４２，４３，・・・それぞれに対して実行したＱＴＳＡの計算結果は基底状態になる場合と励起状態になる場合とがあり得る。そこで初期状態として十分な数のスピン情報４１，４２，４３，・・・を生成することで、いずれかの計算結果は基底状態となることが期待できる。

【0202】

また、このような量子試行とノイズの利用は、計算時間についても短縮させることができる。この理由を以下に説明する。
虚数時間発展において、式（２．５）より、励起状態の減衰速度は式（５．２）を満たすことが要求される。

【0203】

【数68】

【0204】

ここでエネルギーＥ_nは横磁場項を入れ、各々のＧ₁、Ｇ₂に対して定義されたハミルトニアンの式（１．５）に対して定義される。これをＥ_n（Ｇ₁，Ｇ₂）と表記する。このとき、計算の律速になるのは、最も減衰が遅い項である。これは第一励起状態である。従って、計算にかかる時間は凡そ式（５．３）に示す程度である。

【0205】

【数69】

【0206】

ここでＥ₀からＥ₁への変化量ΔＥ₁は式（５．４）で表される。

【0207】

【数70】

【0208】

ここで１次の量子相転移が起こるのであれば、式（５．５）に示す程度のエネルギーギャップが存在する。

【0209】

【数71】

【0210】

従って式（５．６）に示す程度の計算時間を要する。

【0211】

【数72】

【0212】

これは計算時間が実質的に指数関数的に延びることを意味する。このような計算時間は不確定性原理から決まってしまうため、回避するのは難しい。
しかし、実質的に縮退が起きていると決めてしまえば、エネルギーの打切り誤差の範囲で計算量の増大を抑えることが可能である。その代わり、そのエネルギーの打切り誤差の範囲内にある全ての量子状態が候補となるが、これは全体の解空間の大きさと比べれば小さくなる。なぜなら、全ての量子状態が縮退することは実質的にないからである。

【0213】

以上説明したように、最適化装置１００は、スピンの初期状態にノイズを付加してＱＴＳＡ法による最小値求解問題の解を計算することで、効率的に解を求めることができる。
＜ノイズ付加を伴うＱＴＳＡ法の最適化装置への実装＞
次に、Ｑスピンの初期状態にノイズを付加してＱＴＳＡ法により最小値求解問題を解く最適化装置１００の機能について具体的に説明する。

【0214】

図８は、最適化装置の機能を示すブロック図である。最適化装置１００は、記憶部１１０と処理部１２０とを有する。記憶部１１０は、例えばメモリ１０２またはストレージ装置１０３の記憶領域の一部によって実現される。処理部１２０は、例えばプロセッサ１０１に所定のプログラムを実行させることにより実現される。

【0215】

記憶部１１０は、エネルギー情報、スピン情報、横磁場情報、問題設定情報、ハミルトニアン情報を記憶する。
エネルギー情報は、計算されたエネルギーの初期値や、これまで計算されたエネルギーの値のうち、最小値から小さい順に所定数のエネルギーの値を含む。また、エネルギー情報は、所定数のエネルギーの値に対応した各状態変数の値の組合せを含んでいてもよい。スピン情報は、各状態変数の値や、スピン初期化手法の指定情報を含む。横磁場情報は、系に作用する横磁場の強度の情報を含む。問題設定情報は、例えば、求解対象の問題を示す情報、使用する最適化方法（ＳＡ法、ＱＴＳＡ法など）、最適化方法の実行に用いる各種パラメータ、計算終了条件などを含む。ハミルトニアン情報は、例えば、式（１．１）に示したエネルギー関数の重み係数（Ｗ_ij）、バイアス係数（ｂ_i）、定数（Ｃ）などを含む。

【0216】

処理部１２０は、制御部１２１、設定読込部１２２、スピン初期化部１２３、ハミルトニアン生成部１２４、横磁場計算部１２５、虚数時間発展制御部１２６、ノイズ付加部１２７、および古典状態化部１２８を有する。

【0217】

制御部１２１は、処理部１２０の各部を制御して最小値求解処理を行う。
設定読込部１２２は、記憶部１１０から上記の各種情報を、制御部１２１が理解可能な形式で読み込む。

【0218】

スピン初期化部１２３は、イジングモデルに含まれる複数のスピンの状態を初期化する。スピン初期化部１２３は、初期化されたスピンの状態を示す状態変数の値を含むスピン情報を、記憶部１１０に格納する。

【0219】

ハミルトニアン生成部１２４は、問題設定情報に基づいて、求解対象の問題を求解するためのイジングモデルのハミルトニアンを示すハミルトニアン情報を生成する。ハミルトニアン生成部１２４は、生成したハミルトニアンを記憶部１１０に格納する。

【0220】

横磁場計算部１２５は、実数時間の発展に伴って更新される横磁場を計算する。
虚数時間発展制御部１２６は、虚数時間の発展を制御する。例えば虚数時間発展制御部１２６は、横磁場を固定した状態でスピンの状態が計算されるごとに、虚数時間の時間刻み幅Δτずつ虚数時間を進める。

【0221】

ノイズ付加部１２７は、スピン初期化部１２３で初期化されたスピン状態にノイズを付加する。例えばノイズ付加部１２７は、複数のスピンのうちの所定数のスピンをランダムに選択する。ノイズ付加部１２７は、選択したスピンの確率振幅に対してノイズを付加する。ノイズ付加部１２７は、スピンごとのノイズの量を、例えば乱数を用いて個別に決定する。

【0222】

古典状態化部１２８は、量子試行から得られた解の集合から古典イジングモデルに対応する状態を求め、古典ＭＣＭＣ解法の初期状態にする。そして古典状態化部１２８は、古典ＭＣＭＣにより基底状態を求める。

【0223】

なお、図８に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。また、図８に示した各要素の機能は、例えば、その要素に対応するプログラムモジュールをコンピュータに実行させることで実現することができる。

【0224】

次に最適化装置１００による、ＱＴＳＡ法を用いた最小値求解問題の解探索手順について説明する。
図９は、解探索手順の一例を示すフローチャートである。以下、図９に示す処理をステップ番号に沿って説明する。

【0225】

［ステップＳ１０１］制御部１２１は、記憶部からハミルトニアン情報の読み込み処理を行う。この処理の詳細は後述する（図１０参照）。
［ステップＳ１０２］スピン初期化部１２３は、スピン情報を初期化する。例えばスピン初期化部１２３は、ユーザから指定されたスピン初期化手法を示す指定情報を記憶部１１０から読み込む。スピン初期化手法には、例えば指定モード、結合モード、反結合モード、ランダムモードがある。指定モードは、ユーザから指定された通りのスピンの状態に設定する初期化手法である。結合モードは、各スピンにおける｜１〉の状態と｜０〉の状態との係数の符号を同じにする初期化手法である。反結合モードは、各スピンにおける｜１〉の状態と｜０〉の状態との係数の符号を逆にする初期化手法である。スピン初期化部１２３は、ユーザから指定モード、結合モード、反結合モードのいずれかのスピン初期化手法の指定がある場合、指定された初期化手法でスピン情報を初期化する。ランダムモードは、すべてのスピンについて、確率的に上向きまたは下向きに設定する初期化手法である。ランダムモードの場合、スピン初期化部１２３は、上向きの確率と下向きの確率とを５０％ずつにすることで、上向きのスピンと下向きのスピンとを半分ずつとすることができる。スピン情報初期化処理の詳細は後述する（図１１参照）。

【0226】

［ステップＳ１０３］制御部１２１は、ステップＳ１０４～Ｓ１０５の処理を、終了条件が満たされるまで繰り返す。例えば制御部１２１は、予め設定された繰り返しの計算回数をｔｍａｘに設定する。そして制御部１２１は、繰返し回数を示す変数ｔｓｔｅｐを１から順に、ステップＳ１０４～Ｓ１０５の処理を１回実行するごとに１ずつカウントアップする。制御部１２１は、ｔｓｔｅｐの値がｔｍａｘ以下の間、ステップＳ１０４～Ｓ１０５の処理を繰り返し実行する。

【0227】

［ステップＳ１０４］横磁場計算部１２５は、横磁場項の係数Ｇ₂の値を決定する。例えば横磁場計算部１２５は、Ｇ₂の初期値を「１」とし、ステップＳ１０４～Ｓ１０５の処理を繰り返すごとに、Ｇ₂の値を所定量ずつ低下させる。横磁場計算部１２５は、ｔｓｔｅｐの値がｔｍａｘとなったときには、Ｇ₂の値が「０」になるようにする。

【0228】

［ステップＳ１０５］虚数時間発展制御部１２６、ノイズ付加部１２７、および古典状態化部１２８は互いに連係して量子試行計算を行う。量子試行計算により、現在の横磁場におけるハミルトニアンの基底状態が求められる。量子試行計算の詳細は後述する（図１２参照）。

【0229】

［ステップＳ１０６］制御部１２１は、繰り返し処理の終了条件が満たされた場合、処理をステップＳ１０７に進める。例えば制御部１２１は、ステップＳ１０５の処理の終了時点でｔｓｔｅｐの値がｔｍａｘと等しい場合、終了条件が満たされたと判断する。

【0230】

［ステップＳ１０７］制御部１２１は、終了状態を出力する。例えば制御部１２１は、横磁場の係数Ｇ₂の値が「０」になったときのハミルトニアンの基底状態（最適値）および基底状態における状態変数の組み合わせを出力する。

【0231】

このようにして、虚数時間発展処理を繰り返しながら徐々に横磁場の強さを弱めていくことで、系の最適値を得ることができる状態変数の組み合わせを求めることができる。得られた状態変数の組み合わせが、最小値求解問題の解である。

【0232】

次に、ハミルトニアン情報読み込み処理について詳細に説明する。
図１０は、ハミルトニアン情報読み込み処理の手順の一例を示すフローチャートである。以下、図１０に示す処理をステップ番号に沿って説明する。

【0233】

［ステップＳ２０１］制御部１２１は、求解対象のハミルトニアンの係数Ｗ_ji，ｂ_iおよび定数Ｃを、記憶部１１０から読み込む。この際、制御部１２１は、ＱＴＳＡ法の計算に使用する各種パラメータも記憶部１１０から読み込む。例えば制御部１２１は、虚数時間刻み幅Δτが予め指定されている場合、Δτの値を読み込む。なお、虚数時間の刻み幅Δτとして可変時間刻みを採用することも可能であり、その場合にはΔτを所定の計算で求めることとなり、虚数時間刻み幅Δτの読み込みは不要である。

【0234】

［ステップＳ２０２］スピン初期化部１２３は、スピン初期化手法の指定情報を読み込む。
［ステップＳ２０３］スピン初期化部１２３は、スピン初期化手法が指定モードか否かを判断する。スピン初期化部１２３は、指定モードであれば、処理をステップＳ２０４に進める。またスピン初期化部１２３は、指定モードでなければ、ハミルトニアン情報読み込み処理を終了する。

【0235】

［ステップＳ２０４］スピン初期化部１２３は、ユーザによって指定されたスピンの値を記憶部１１０から読み込む。その後、スピン初期化部１２３は、ハミルトニアン情報読み込み処理を終了する。

【0236】

次にスピン情報初期化処理について詳細に説明する。
図１１は、スピン情報初期化処理の手順の一例を示すフローチャートである。以下、図１１に示す処理をステップ番号に沿って説明する。

【0237】

［ステップＳ３０１］スピン初期化部１２３は、スピン初期化手法が指定モードか否かを判断する。スピン初期化部１２３は、指定モードであれば、処理をステップＳ３０２に進める。またスピン初期化部１２３は、指定モードでなければ、処理をステップＳ３０３に進める。

【0238】

［ステップＳ３０２］スピン初期化部１２３は、ユーザによって指定された値にスピンの状態を初期化する。例えばスピン初期化部１２３は、記憶部１１０から読み込んだスピンの状態の情報に基づいて、スピンの状態を初期化する。スピン初期化部１２３は、その後、スピン情報初期化処理を終了する。

【0239】

［ステップＳ３０３］スピン初期化部１２３は、スピン初期化手法が結合モードか否かを判断する。スピン初期化部１２３は、結合モードであれば、処理をステップＳ３０４に進める。またスピン初期化部１２３は、結合モードでなければ、処理をステップＳ３０５に進める。

【0240】

［ステップＳ３０４］スピン初期化部１２３は、結合モードとなるようにスピンの状態を初期化する。例えばスピン初期化部１２３は、すべてのスピンの状態を「｜φ_i〉＝１／（２^1/2）（｜１〉＋｜０〉）」と初期化する。スピン初期化部１２３は、その後、スピン情報初期化処理を終了する。

【0241】

［ステップＳ３０５］スピン初期化部１２３は、スピン初期化手法が反結合モードか否かを判断する。スピン初期化部１２３は、反結合モードであれば、処理をステップＳ３０６に進める。またスピン初期化部１２３は、反結合モードでなければ、処理をステップＳ３０７に進める。

【0242】

［ステップＳ３０６］スピン初期化部１２３は、反結合モードとなるようにスピンの状態を初期化する。例えばスピン初期化部１２３は、すべてのスピンの状態を「｜φ_i〉＝１／（２^1/2）（｜１〉－｜０〉）」と初期化する。スピン初期化部１２３は、その後、スピン情報初期化処理を終了する。

【0243】

［ステップＳ３０７］スピン初期化部１２３は、ランダムモードとなるようにスピンの状態を初期化する。例えばスピン初期化部１２３は、すべてのスピンの状態のうちの確率５０％で選択されたスピンを「｜φ_i〉＝｜０〉」、残りを「｜φ_i〉＝｜０〉」と初期化する。スピン初期化部１２３は、その後、スピン情報初期化処理を終了する。

【0244】

このようにしてスピン情報を初期化することができる。初期化されたスピン情報に基づいて設定される量子ビット（スピン）ごとの波動関数の一部には、量子試行計算の初期段階でノイズが付加される。以下、ノイズ付加を伴う量子試行計算処理の手順について詳細に説明する。

【0245】

図１２は、量子試行計算処理の手順の一例を示すフローチャートである。以下、図１２に示す処理をステップ番号に沿って説明する。
［ステップＳ４０１］虚数時間発展制御部１２６は、量子試行回数Ｎ_qなどの計算パラメータを設定する。例えば虚数時間発展制御部１２６は、予め指定されている計算パラメータの値を対応する変数に設定する。

【0246】

［ステップＳ４０２］虚数時間発展制御部１２６は、ステップＳ４０３～Ｓ４０９の処理をＮ_q回繰り返す。
［ステップＳ４０３］虚数時間発展制御部１２６は、スピンごとの波動関数の初期状態を設定する。例えば虚数時間発展制御部１２６は、スピン情報初期化処理によって初期化された値に応じた式（４．１）の波動関数をスピンごとに生成する。

【0247】

［ステップＳ４０４］ノイズ付加部１２７は、波動関数へのノイズ付加処理を行う。ノイズ付加処理の詳細は後述する（図１３参照）。
［ステップＳ４０５］虚数時間発展制御部１２６は、虚数時間発展を用いたＱＡを実行する。虚数時間発展を用いたＱＡ処理の詳細は後述する（図１４参照）。

【0248】

［ステップＳ４０６］古典状態化部１２８は、現在のスピンの状態の対応する古典的状態を求める。古典的状態を求める処理では、古典的状態が存在しない場合、処理結果として準安定状態「metastable state」が出力される。古典的状態を求める処理の詳細は後述する（図１５参照）。

【0249】

［ステップＳ４０７］古典状態化部１２８は、現在のスピンの状態に対応する古典的状態が存在するか否かを判断する。例えば古典状態化部１２８は、古典的状態を求める処理結果が「metastable state」となった場合、古典的状態が存在しないと判断し、それ以外の場合は古典的状態が存在すると判断する。古典状態化部１２８は、対応する古典的状態が存在する場合、処理をステップＳ４０８に進める。また古典状態化部１２８は、対応する古典的状態が存在しない場合、処理をステップＳ４０９に進める。

【0250】

［ステップＳ４０８］古典状態化部１２８は、ステップＳ４０６で得られた古典的状態を初期状態として、古典ＭＣＭＣにより基底状態を求める。例えば古典状態化部１２８は、古典状態のイジングモデルのどのスピンを反転対象とするか、およびそのスピンの反転を受け入れるか否かについて乱数を用いて確率的に決定する。古典状態化部１２８は、スピンの反転を受け入れると決定した場合、そのスピンを反転した状態にイジングモデルの状態を遷移させる。古典状態化部１２８は、このような状態遷移を繰り返してエネルギーが最小となる状態を求める。そして古典状態化部１２８は、古典ＭＣＭＣで得られたエネルギー値とスピン状態を現在の処理ループでの処理結果として保存する。古典状態化部１２８は、その後、処理をステップＳ４１０に進める。なお、古典ＭＣＭＣによる求解は非特許文献１に記載されているような装置と連携させて求めることもできる。

【0251】

［ステップＳ４０９］古典状態化部１２８は、対応する古典状態は存在しないとの処理結果を出力する。
［ステップＳ４１０］虚数時間発展制御部１２６は、ステップＳ４０３～Ｓ４０９の処理のＮ_q回の繰り返しが完了した場合、量子試行計算処理を終了する。

【0252】

［ステップＳ４１１］虚数時間発展制御部１２６は、ステップＳ４０８で保存したエネルギー値のうちの最小のエネルギー値と、最小のエネルギー値を得たときのスピン状態を、量子試行計算の結果として出力する。

【0253】

このようにノイズが付加された波動関数を初期状態とする虚数時間発展を用いたＱＡが複数回実行され、エネルギー値が最小となるスピン状態が得られる。波動関数へのノイズの付加は乱数を用いて確率的に行われる。

【0254】

図１３は、波動関数へのノイズ付加処理の手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。
［ステップＳ５０１］ノイズ付加部１２７は、ノイズを載せるスピン数Ｍを設定する。例えばノイズ付加部１２７は、予め指定された値を、ノイズを載せるスピン数Ｍを示す変数に設定する。

【0255】

［ステップＳ５０２］ノイズ付加部１２７は、ノイズの最大の大きさｚ_max［％］を設定する。例えばノイズ付加部１２７は、予め指定された値を、ノイズの大きさｚを示す変数に設定する。

【0256】

［ステップＳ５０３］ノイズ付加部１２７は、ステップＳ５０４～Ｓ５０７の処理をＭ回繰り返す。
［ステップＳ５０４］ノイズ付加部１２７は、乱数を用いて１からＮまでのうちの１つの整数ｊを発生させる。

【0257】

［ステップＳ５０５］ノイズ付加部１２７は、最後に生成した乱数が、それ以前に生成した乱数と重複するか否かを判断する。ノイズ付加部１２７は、乱数の重複がある場合、処理をステップＳ５０４に進める。またノイズ付加部１２７は、乱数の重複がなければ、処理をステップＳ５０６に進める。

【0258】

［ステップＳ５０６］ノイズ付加部１２７は、乱数を用いて０からｚ_max［％］までの実数値を発生させる。そしてノイズ付加部１２７は、発生した実数値をｚ［％］とする。
［ステップＳ５０７］ノイズ付加部１２７は、ｊ番目のスピンの波動関数における確率振幅にｚ［％］のノイズを加える。

【0259】

［ステップＳ５０８］ノイズ付加部１２７は、ステップＳ５０４～Ｓ５０７の処理のＭ回の繰り返しが完了した場合、波動関数へのノイズの付加処理を終了する。
このようにしてＭ個のスピンについてランダムな大きさのノイズが付加された波動関数について、虚数時間発展処理を用いたＱＡが実行される。なおノイズの大きさを一律にしてもよい。その場合、ノイズ付加部１２７は、例えばステップＳ５０６における乱数による実数値の発生処理を行わずに、ｚ_max［％］をｚ［％］とする。

【0260】

図１４は、虚数時間発展を用いたＱＡ処理の一例を示すフローチャートである。以下、図１４に示す処理をステップ番号に沿って説明する。なお、以下の説明では、上記の数式に示されるスピンの番号を示す記号をｉからｋに置き換えている。

【0261】

［ステップＳ６０１］虚数時間発展制御部１２６は、式（１．１９）に基づいて、縦磁場演算子の期待値ｘ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。
［ステップＳ６０２］虚数時間発展制御部１２６は、式（１．２０）に基づいて、横磁場演算子の期待値ｙ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。

【0262】

［ステップＳ６０３］虚数時間発展制御部１２６は、式（２．２３）と式（２．２４）に基づいて、縦磁場平均場ｈ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。

【0263】

［ステップＳ６０４］虚数時間発展制御部１２６は、式（２．３５）に基づいて、横磁場エネルギーＥ_OLを計算する。
［ステップＳ６０５］虚数時間発展制御部１２６は、ステップＳ６０６～Ｓ６１３の処理を、終了条件が満たされるまで繰り返す。例えば虚数時間発展制御部１２６は、予め設定された繰り返しの計算回数をｎｉｔｌに設定する。そして制御部１２１は、虚数時間発展回数を示す変数ｉを１から順に、ステップＳ６０６～Ｓ６０９の処理を１回実行するごとに１ずつカウントアップする。制御部１２１は、ｉの値がｎｉｔｌ以下の間、ステップＳ６０６～Ｓ６０９の処理を繰り返し実行する。

【0264】

［ステップＳ６０６］虚数時間発展制御部１２６は、式（２．１６）に基づいて、次の時刻（Δτだけ虚数時間を進行させた時刻）のＣ_kα（τ₀＋Δτ），Ｃ_kβ（τ₀＋Δτ）を計算する。

【0265】

［ステップＳ６０７］虚数時間発展制御部１２６は、以下の式（６．１）に基づいて、規格化条件を満たすように再規格化を行う。

【0266】

【数73】

【0267】

γ_mは規格化定数である。再規格化前の波動関数の係数Ｃ_mα（ｔ＋Δｔ）、Ｃ_mβ（ｔ＋Δｔ）それぞれに規格化定数γ_mを乗算した結果が、再規格化後の波動関数の係数となる。

【0268】

［ステップＳ６０８］虚数時間発展制御部１２６は、再規格化した係数を用いた式（１．１９）に基づいて、縦磁場演算子の期待値ｘ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。

【0269】

［ステップＳ６０９］虚数時間発展制御部１２６は、再規格化した係数を用いた式（１．２０）に基づいて、横磁場演算子の期待値ｙ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。

【0270】

［ステップＳ６１０］虚数時間発展制御部１２６は、式（２．２３）と式（２．２４）に基づいて、縦磁場平均場ｈ_kをすべてのスピンについて計算する（ｋ＝１，２，・・・，Ｎ）。

【0271】

［ステップＳ６１１］虚数時間発展制御部１２６は、式（２．３５）に基づいて、横磁場エネルギーＥ_OLを計算する。
［ステップＳ６１２］虚数時間発展制御部１２６は、式（２．３４）に基づく全エネルギーＥ（Ｇ₁，Ｇ₂）と、前の時刻との全エネルギーの差ΔＥを計算する。

【0272】

［ステップＳ６１３］虚数時間発展制御部１２６は、前の時刻との全エネルギーの差ΔＥが計算の打ち切り精度ε以下か否かを判断する。虚数時間発展制御部１２６は、ΔＥがε以下の場合、虚数時間発展処理を終了する。虚数時間発展制御部１２６は、ΔＥがεより大きければ処理をステップＳ６１４に進める。

【0273】

［ステップＳ６１４］虚数時間発展制御部１２６は、繰り返し処理の終了条件が満たされた場合、虚数時間発展処理を終了する。例えば虚数時間発展制御部１２６は、ステップＳ６１３の処理の終了時点でｉの値がｎｉｔｌと等しい場合、終了条件が満たされたと判断する。

【0274】

量子試行計算が終了すると、対応する古典的状態を求める処理が実行される。
図１５は、対応する古典的状態を求める処理の手順の一例を示すフローチャートである。以下、図１５に示す処理をステップ番号に沿って説明する。

【0275】

［ステップＳ７０１］古典状態化部１２８は、閾値εを設定する。例えば古典状態化部１２８は、予め指定された値を、閾値εを示す変数に設定する。
［ステップＳ７０２］古典状態化部１２８は、変数ｉを１からカウントアップしながらステップＳ７０３～Ｓ７０９の処理をＮ（スピン数）回繰り返す。

【0276】

［ステップＳ７０３］古典状態化部１２８は、確率振幅を示す値（式（４．１）参照）であるＣ_iαとＣ_iβとの絶対値の大小関係を比較する。古典状態化部１２８は、Ｃ_iαの方が大きければ処理をステップＳ７０４に進める。古典状態化部１２８は、｜Ｃ_iα｜が｜Ｃ_iβ｜以下であれば処理をステップＳ７０７に進める。

【0277】

［ステップＳ７０４］古典状態化部１２８は、｜Ｃ_iα｜が１－εより大きいか否かを判断する。古典状態化部１２８は、｜Ｃ_iα｜が１－εより大きければ、処理をステップＳ７０６に進める。また古典状態化部１２８は、｜Ｃ_iα｜が１－ε以下であれば、処理をステップＳ７０５に進める。

【0278】

［ステップＳ７０５］古典状態化部１２８は、準安定状態「metastable state」である旨を出力し、対応する古典的状態を求める処理を終了する。
［ステップＳ７０６］古典状態化部１２８は、ｉ番目のスピンの状態を１にする。その後、古典状態化部１２８は処理をステップＳ７１０に進める。

【0279】

［ステップＳ７０７］古典状態化部１２８は、｜Ｃ_iβ｜が１－εより大きいか否かを判断する。古典状態化部１２８は、｜Ｃ_iβ｜が１－εより大きければ、処理をステップＳ７０９に進める。また古典状態化部１２８は、｜Ｃ_iβ｜が１－ε以下であれば、処理をステップＳ７０８に進める。

【0280】

［ステップＳ７０８］古典状態化部１２８は、準安定状態「metastable state」である旨を出力し、対応する古典的状態を求める処理を終了する。
［ステップＳ７０９］古典状態化部１２８は、ｉ番目のスピンの状態を０にする。

【0281】

［ステップＳ７１０］古典状態化部１２８は、ステップＳ７０３～Ｓ７０９の処理のＮ回の繰り返しが完了したら、対応する古典状態を求める処理を終了する。
古典状態を求めることで、量子試行だけでは基底状態が求まらない場合であっても、波動方程式の基底状態を求めることができる。その結果、虚数時間発展処理を用いたことにより基底状態が得られなくなる可能性を低減することができる。

【0282】

＜組み合わせ最適化問題の計算例＞
以下、実際の組み合わせ最適化問題の実際の計算例について説明する。なお以下の計算例は、スタンドアロンのノイマン型コンピュータで計算したものである。計算は、並列化を行わず、ＧＰＵなどのアクセラレータも使用せずに、実用的な時間内で実施されている。

【0283】

＜＜初期状態にノイズを付加した量子試行による崩壊現象の実例＞＞
最適化装置１００により２００回の量子試行を行い、前述した４都市ＴＳＰを解いた場合の計算結果について説明する。最適化装置１００に与えられた計算の条件は以下の通りである。

【0284】

最適化装置１００は、任意のｋ（ｋ＝１，２，・・・）スピンを選んでノイズを付加している。また最適化装置１００は、Ｇ₁=１．０と固定し、Ｇ₂を１から０へ変化させて問題を解いている。Ｇ₂の分点数はｎ₀＝１２８である。Ｇ₂（ｎ）は以下の式（７．１）で表される。

【0285】

【数74】

【0286】

最適化装置１００は、ＱＴＳＡ（虚数時間発展）の打切り誤差はδＥ＝１０^-8としている。また最適化装置１００は、与えられたＧ₂（ｎ）においてＱＴＳＡ計算がエネルギーの収束の打切り誤差以下になるまで継続する。最適化装置１００は、エネルギーが収束した時点でｎの値を１つ増やし、再度ＱＴＳＡ計算を繰り返す。最終的にＧ₂=０のときの値が得られている。

【0287】

最適化装置１００は、波動関数の初期値には１６個の量子ビット（スピン）の中から任意の２つの量子ビットをランダムで選択し、選択されたビットの確率振幅Ｃ_mα（０）に最大１％のノイズをランダムに入れて計算する。このような量子試行計算を２００回ほど繰り返した結果を図１６に示す。

【0288】

図１６は、４都市ＴＳＰにおける２００回の量子試行結果の一例を示す図である。図１６において、左側のグラフ５１は、横磁場の強さに応じたエネルギーの遷移を示している。グラフ５１の横軸が横磁場の強さであり、縦軸がエネルギーである。右側のグラフ５２は、虚数時間（Ｉｍａｇｉｎａｒｙ Ｔｉｍｅ）の進行に応じたエネルギーの遷移を示している。グラフ５２の横軸が虚数時間であり、縦軸がエネルギーである。

【0289】

グラフ５１，５２には、２００回の量子試行それぞれに応じた実線が重ねて描かれている。グラフ５１を参照すると、すべての試行において、４都市ＴＳＰの解であるエネルギー状態７．４１４，７．８１２，１０．３９８のいずれかに到達していることが分かる。基底状態に１００％到達するわけではないが、有限の確率で到達していることが分かる。なお、解に到達するまでの過程は試行により異なるが、多数の試行においてエネルギーが高い状態を経由していることが分かる。

【0290】

エネルギーが高い量子状態の特徴はグラフ５２に現れている。特徴的であるのは、エネルギーの高い状態が有限の寿命を持ってＴＳＰの解の状態に崩壊（ｃｏｌｌａｐｓｅ）していることである。これは準安定状態と呼ばれるものである。

【0291】

虚数時間発展を続けると、いきなり最終状態に行き着く場合もあるが、準安定状態を経由する場合もある。そして、この準安定状態は一般に１と０のみで記述される古典的な状態とは限らない。

【0292】

図１７は、スピンの波動関数の時間変化の一例を示す図である。左のグラフ５３は、準安定状態を経由せず基底状態に到達した場合の１６個のスピンの波動関数の時間変化を示している。右のグラフ５４は、準安定状態を経由して基底状態に到達した場合の１６個のスピンの波動関数の時間変化を示している。グラフ５３，５４は横軸が虚数時間であり縦軸が波動関数の絶対値の２乗である。

【0293】

準安定状態を経由しない場合、虚数時間発展の初期ですべてのスピンが反転するのに対し、準安定状態を経由する場合は、虚数時間発展の初期で１６個のスピンのうち１２個は反転するが、残りの４つのスピンの役割が決まらず、中間的な状態を取っている。この状態は対応する古典的な状態は存在しない。量子力学特有の状態である。そして、虚数時間τ≒６程度で４つのスピン中の２つペアが０と１にそれぞれ分かれて基底状態に到達する。これは一種の分岐現象（ｂｉｆｕｒｃａｔｉｏｎ）である。

【0294】

なお、系の量子ビット数が大きくなると準安定状態の寿命が延びる傾向にある。これは深いポテンシャル障壁の存在を示唆するが、寿命が長く安定であるため、熱散逸をもってしてもエネルギーが低い状態に崩壊せず計算が終了する場合もある。この場合、対応する古典的状態は存在しないため、対応する古典イジングモデルの解を求めることができない。そのため最適化装置１００は、古典イジングモデルの解を抜き出す際に、準安定状態で終わってしまった解を準安定状態（metastable）として答えの候補から外す。しかし準安定状態は古典論と対応が付かないだけであり、その状態が無意味であるというわけではない。本質的に量子論の方が古典論の拡張になっているため生じるのである。

【0295】

＜＜ノイズ付加による対称性の破壊の実例＞＞
次に、スピンの波動関数の初期値へノイズを与えることで意図的に波動関数の対称性の破壊を誘導可能であることの実例を示す。計算対象の問題は８都市のＴＳＰである。これは６４量子ビットの系である。

【0296】

図１８は、８都市のＴＳＰの一例を示す図である。図１８には、各都市の識別番号に対応付けて、その都市の位置（ｘ座標とｙ座標）が示されている。この問題の最適解の１つでは、都市の訪問順序は都市１→２→３→４→５→６→７→８である。このときの最適値は－２．５０２３である。

【0297】

図１９は、８都市のＴＳＰの量子試行結果の一例を示す図である。図１９に示すグラフ５５は、横軸が横磁場の強さであり、縦軸がエネルギーである。グラフ５５に示される各線は、異なるノイズが付加された複数の初期状態それぞれに基づく量子試行における横磁場の強さに応じたエネルギーの遷移を示す。

【0298】

初期状態にノイズを入れたため、式（３．２）に従い散乱され終状態に達するが、初期状態が異なるため異なる終状態に達している。このため、基底状態を含む様々な量子状態が得られることが分かる。このように小さいノイズであっても初期状態に加わると、量子アニーリングの結果として必ずしも基底状態が得られないことが分かる。

【0299】

これは１００％の確率で基底状態を得たい場合には望ましくない挙動である。しかし、ノイズの効果を逆手に取れば、ノイズを意図的に加えれば基底状態を含む量子状態を計算できることを意味している。そして、基底状態も１００％ではないが有限の確率で計算できる。

【0300】

次にノイズ付加により波動関数の対称性が破壊される様子を示す。
図２０は、８都市ＴＳＰにおける波動関数の時間発展の一例を示す図である。図２０には、基底状態と第一励起状態について６４個の量子ビットの波動関数の時間発展の様子が示されている。左側のグラフ５６が基底状態の時間発展を示し、右側のグラフ５７が第１励起状態の時間発展を示す。グラフ５６，５７では、６４個の｜０＞の時間発展が実線で表され、６４個の｜１＞の時間発展が点線で表されている。グラフ５６，５７は、横軸が虚数時間であり、縦軸が波動関数の絶対値の２乗である。

【0301】

基底状態と第一励起状態の比較から分かるように、ノイズを付加しても、系の時間発展の初期段階ではあまり違いがない。しかし、時間が経つにつれノイズによる影響が増幅されて行く。そして、結果としてどのスピンが反転するかという系全体の振る舞いにも影響を与えてしまうのである。このため、波動関数の初期値にわずかにノイズを乗せるだけで、もはや決定論的には基底状態に到達しないのである。

【0302】

次にスピン反転の様子について説明する。時間発展の初期段階ではすべてのスピンが協調的に足並みをそろえて反転していく。しかし、どれか一つのスピンが他のスピンより高速に反転を開始すると、系全体の対称性が崩れる。そして、反転したスピンの影響により別のスピンが反転を開始するという具合に系全体のスピン状態が決定されていく。

【0303】

また、一度対称性が壊れた後に、再度ノイズを乗せても無駄である。スピンを完全に反転させるような処置をすれば別であるが、反転を開始したあとは、多少のノイズを乗せても系はびくともしない。これは物理的には重要である。なぜならば、波動関数の対称性が計算の安定、不安定、量子アニーリングの成功不成功に影響を与えるからである。スピン情報の初期化直後は、例えばすべてのスピンの波動関数が同じ状態になるため、高い対称性があるが、これは容易に壊れる対称性である。半面、一度対称性が壊れてしまうと、多少のノイズを加えても波動関数の対称性を壊してまで低エネルギー状態に移れないのである。そのため、量子アニーリングにおいて１００％の確率で基底状態を求めたいならば、初期状態から終状態まで基底状態の対称性が選択されるまではノイズを乗せてはいけないのである。そのため、断熱的に遷移をさせるなどの系の遷移に対して厳しい制約が課される。

【0304】

一方で、熱的散逸を起こすことを前提として、ノイズを積極的に利用するならば基底状態に到達する確率は１００％にはならないが、代わりに励起状態まで含めた状態を求めることができる。

【0305】

＜＜量子試行と古典ＭＣＭＣの試行の効率の比較＞＞
次にノイズを付加した量子試行で得られた結果と、同じ問題を古典イジングモデルとして扱い、レプリカ交換法の下での古典ＭＣＭＣの計算結果とを比較する。これは量子試行によるサンプリング計算と、古典ＭＣＭＣベースのサンプリング計算の効率を比較するために実施したものである。最適化装置１００に与えられた計算の条件は以下の通りである。

【0306】

古典イジングモデルにおける数値解法はＭＣＭＣベースの解法である。解法には２種類用いている。１つ目はＮ個あるスピンのうち１個をランダムに選びスピンを反転させる試行である。２つ目は１のスピンの数を保存するようにＮ個のスピンの中から１のスピンをランダムに１つ選び０にし、同時に０のスピンをランダムに１つ選択して１にする操作を同時に行うペアフリップである。

【0307】

ペアフリップを行ったのはＴＳＰの場合１のスピンの数が都市数になると決まっているためである。ペアフリップとしたため、答えの数が保存される試行になっており、解の探索空間がシングルフリップよりも大幅に減る。反転の基準はメトロポリス法が適用されている。さらに、単純なＭＣＭＣでは遅すぎるため、レプリカ交換法を用いて古典ＭＣＭＣ計算が行われている。レプリカ数は１０個に固定している。１レプリカあたり１０億回の試行を行い、レプリカ交換頻度は１００回に１度としている。温度については手動で最適化を行っている。

【0308】

図２１は、８都市ＴＳＰにおける量子試行と古典イジングモデルのＭＣＭＣとの計算結果の比較例を示す図である。図２１のグラフ５８には、８都市ＴＳＰにおける量子試行１０００回のうちエネルギーが低い方から上位１００までの結果が丸印で示されている。また、８都市ＴＳＰにおける古典イジングモデルのＭＣＭＣにより得られたエネルギーのうちの低い方から上位１００までの結果が四角の印で示されている。グラフ５８は、横軸がランキング（エネルギーが低い方から何番目か）であり、縦軸がエネルギーである。

【0309】

基底状態については量子論、古典論の場合両方とも得られている。また、励起状態についても下から１１番目（量子論で順位１００位まで）までは一致している。基底状態については８重に縮退しているが、８個以上の基底状態が得られているのはまったく同じ量子状態が求まっているからである。これより、本手法により基底状態を含めた励起状態も求めることが可能であることが分かる。

【0310】

図２２は、試行回数に応じたエネルギーの遷移例を示す図である。図２２のグラフ５９は、試行回数の関数としてエネルギー値が更新される様子を示している。グラフ６０は、グラフ５９のエネルギーが０～５の範囲のみを拡大したものである。グラフ５９，６０は、横軸が試行回数であり、縦軸がエネルギーである。

【0311】

量子試行（丸印を結ぶ折れ線）についてはこの図２２の例では７５回目で基底状態の値が得られている。ただし、これは乱数の与え方によって異なり、乱数シードを変えた場合３回目で基底状態に到達する例も得られている。一方で古典イジングをＭＣＭＣベースで解いた場合（三角の印を結ぶ折れ線）、４７４９２回目の試行で解にたどり着いている。レプリカ数を考慮すると４７４９２０回目の試行で基底状態にたどり着いたことになる。

【0312】

しかし、この４７４９２０回という数値の絶対値についてはあまり意味がないことを留意すべきである。なぜなら、ＭＣＭＣベースのレプリカ交換は最適化によって１０倍から１００倍加速できるからである。ＭＣＭＣベースの解法で決まっているのは最小値の更新の仕方が初期のバーンイン（初期緩和）が終わった後、ほぼ対数の逆数に比例するように非常にゆっくり更新していくことである。そのため、スケーリング特性が非常に重要である。

【0313】

８都市ＴＳＰの場合、このスケーリング特性が顕著にはなっていないが、都市数が増えるとスケーリング特性の違いが顕著になっていく。これは後で都市数を増やした場合の計算結果で示す。

【0314】

量子論の結果として重要なことは、ノイズと熱散逸の効果により脱励起を起こすことである。脱励起はトンネル効果で高エネルギー状態から低エネルギー状態へ遷移していると考えてもよい。このため、終状態も基底状態に比較的近い状態が得られる。

【0315】

＜＜ulysses16の場合におけるスケーリング評価＞＞
次に古典ＭＣＭＣベースの解法と比較した場合の解の更新速度のスケーリングを比較する。問題にはＴＳＰＬＩＢに収録されているulysses16を用いる。これは最適解が６８５９と知られている問題である。そしてこの問題は８都市のＴＳＰと異なり巡回経路はかなり多くなるため、古典論と量子論の違いが明確に出せるほどには複雑な問題である。

【0316】

ここで８都市の場合と同様にして、古典イジングモデルの解法はレプリカ交換法で１０レプリカを使い、１レプリカあたり１０億回計算を行った計算例について説明する。この計算例では、レプリカ交換頻度は１００回に一度である。ＱＴＳＡはＧ₁=０．００１、Ｇ₂は１から０まで式（６．１）に従って減らされている。量子試行の回数は１０００回である。ノイズは１粒子励起的に入れられている。

【0317】

図２３は、ulysses16におけるＱＴＳＡと古典ＭＣＭＣで得られた解の比較例を示す図である。左側のグラフ６１は、古典ＭＣＭＣとＱＴＳＡの結果のエネルギーが低い方から上位１００位までを比較した図である。グラフ６１の横軸は順位であり、縦軸はエネルギーである。右側のグラフ６２は、横磁場の強さＧ₂の関数としてエネルギー（求解対象）の変化を図示したものである。グラフ６２の横軸は横磁場の強さであり、縦軸はエネルギーである。

【0318】

図２３を参照すると、まず、量子論に基づく計算（グラフ６１では丸印）でも基底状態に到達できていることが分かる。次に古典ＭＣＭＣの計算結果（グラフ６１では四角または△の印）との比較においては、量子論に基づく計算のほうが良い結果が得られている。古典ＭＣＭＣではレプリカごとに１０億回の計算をするこのパラメータセットでは解に到達していない。なおグラフ６１におけるＳＳＦはシングルスピンフリップ（Single Spin Flip）、ＰＳＦはペアスピンフリップ（Pair Spin Flip）を表す。ＳＳＦの方は最小値が７４４２、ＰＳＦの最小値は７１８１である。量子論で得られた最小値は６８５９であり、公表されている最小値と一致している。

【0319】

図２４は、解の改善速度の量子論と古典論との比較例を示す図である。図２４では、解の更新の様子を試行回数の関数として図示している。左側のグラフ６３は、１ステップ目からの解の更新の様子を示している。右側のグラフ６４は、バーンインが終了した後の解の更新の様子を示している。グラフ６３，６４は、横軸が試行回数であり、縦軸がエネルギーである。

【0320】

グラフ６３における最初の１０００ステップはバーンインであるとみなしてよい。バーンインの部分の更新は求解更新速度について本質的ではない。グラフ６３，６４に示すように、解の更新は古典ＭＣＭＣベースの解法（四角または三角の印を結ぶ折れ線）の場合、初期のバーンインが終了した後は更新速度が極端に遅くなり、計算回数に関して対数的にゆっくり更新されている。対して量子論の方（丸印を結ぶ折れ線）はいきなり８，０００付近からデータ計測が始まり、５９７回目で解に到達している。

【0321】

次に、分岐現象が計算を加速するメカニズムについて説明する。この現象についてはスピンの微視的状態の時間変化を見ると説明可能である。
図２５は、典型的な量子試行の時間発展の様子を示す図である。図２５には、ＱＴＳＡ実行時の１つの典型的な量子試行の微視的スピン状態の時間変化が示されている。グラフ６５，６６の横軸が虚数時間であり、縦軸が｜１＞のスピンの観測確率（波動関数の絶対値の２乗）である。グラフ６６は、グラフ６５の観測確率０～０．２の範囲を拡大したものである。

【0322】

グラフ６５からスピン反転は同時には起こらず、特定のスピンが反転して状態が決まった後、連鎖的にスピンが次々に反転していくことが分かる。グラフ６５では虚数時間τ≒１５付近で先に説明した分岐現象が起こる。この現象が起こると、古典イジングモデルでは解候補が２^N個あるものが、１／４の２^N-2になる。なぜならば２つのスピンの古典的状態が確定するためである。もう１つは１つのスピン状態が決まることにより多数の量子ビットが束になって値が古典的状態に確定する場合がある。これがグラフ６６に現れている。１つの量子ビットが１に確定すると多数の量子ビットが束になり、ほぼ同時に０に確定する。この束を構成する量子ビットの数は様々な値をとる。２つであるときもあれば、１０個以上のときもある。重要なことはこの現象が起こると、束を構成する量子ビットの数がｑ個であるなら、その瞬間に解の探索空間は２^N-qになるということである。解空間はその都度指数関数的に減る。例えば、グラフ６６の場合、量子ビットの束が一度に２０個程度一斉に０になっている。２０個程度の量子ビットが｜０＞になると、２²⁰がおよそ１００万分の１であるから、この現象が１回起こるごとに解空間は１００万分の１になっている。この現象が立て続けに短時間で起こるため、解空間は短時間で小さくなる。

【0323】

これは２^Nの広さの解空間を二分木よりも高速な探索をしているに等しい。すなわち、指数関数的に大きい解空間が指数関数的に狭められている。
次に熱散逸によりどのようにエネルギーが落ちていくかについて説明する。

【0324】

図２６は、量子論におけるエネルギー散逸の様子を示す図である。図２６では、ulysses16の最適解である６８５９に到達したときの虚数時間発展による熱散逸の様子が示されている。グラフ６７，６８の横軸が虚数時間であり、縦軸がエネルギーである。グラフ６８は、グラフ６７におけるエネルギーの６８４０～７０００の範囲を拡大したものである。

【0325】

図２６に示す現象が１回ごとの量子試行で必ず起きる。当初横磁場の基底状態に近い状態のときは緩やかにエネルギーが下がっていくが、虚数時間発展が進行すると、ある時間帯において非常に短時間でエネルギーが急激に下がる。この時系の対称性がほぼ決まる。つまり、このときにどのスピンが１になるか、０になるかが決定される。その後、エネルギーがほぼ一定になる時間帯が存在する。図２５のスピン反転図に示されるように、すべてのスピンが一斉に反転するわけではない。一部のスピンの挙動が決まっていない状況になる。このため、ここでは非常にゆっくりエネルギー緩和が起きる。しかし、ある所から最終的に１になるスピンが決定的になると他のスピンも追随して０になる。そして、エネルギーは再度指数関数的に下がり、最終状態である６８５９の状態に行き着く。これは一種の基底選択が行われていると考えると理解しやすい。

【0326】

この傾向は量子ビット数をさらに増やすとより顕著になる。そこでulysses16よりも都市数を増やしたulysses22の場合の例を以下に示す。ulysses22も公開されている問題であり、最適解は７０１３であることが知られている。イジングモデルにすると量子ビット数は４８４Ｑｂｉｔとなる。

【0327】

図２７は、ulysses22におけるランキング情報の一例を示す図である。図２７に示すグラフ６９は、ＱＴＳＡに基づく量子試行２００回で得られたエネルギー値の上位５０位までのエネルギーを図示している。最適化装置１００に与えられた計算の条件は以下の通りである。

【0328】

古典ＭＣＭＣベースの解法におけるレプリカ数は２０個である。レプリカ交換はulysses16と同様にして温度空間をランダムウォークし、十分高い温度まで到達できるようにして計算が実施されている。レプリカごとの計算回数は１０億回である。スピンフリップの仕方はＳＳＦ（グラフ６９中の三角の印）とＰＳＦ（グラフ６９中の四角の印）の両方が実施されている。ペナルティ値は４０００である。

【0329】

グラフ６９中の黒塗りの三角の印は制約条件を満たす解である。ＳＳＦの場合、エネルギーは下がるが、ペナルティ値が４０００であるため、制約違反解を排除できず、得られたエネルギー最小値が７８６１であるが、これは制約違反解である。得られた解の集合の上位５０個の中から制約条件を満たす解だけを探し出したものが、グラフ６９中の黒い三角の印で示されている。制約を満たす解の最小値は８３６５である。従って、最適解までまだエネルギーは１３００以上下がる余地がある。ＰＳＦは初めから制約条件を満たすようにしてスピンフリップをしているため、すべての解が制約を満たす。得られた解の最良値は７８５２である。

【0330】

一方で量子論による結果の最小値は２００回の試行で７０１９である。これは最適解７０１３にはたどり着いていないが、最適解に十分近い値である。
図２８は、ulysses22における解の更新の様子を示す図である。グラフ７０，７１は、横軸が試行回数であり、縦軸がエネルギーである。量子試行の結果が丸印で示され、古典ＭＣＭＣの結果が四角の印で示されている。グラフ７１は、グラフ７０における７０００～１００００の範囲を拡大したものである。

【0331】

自由度が大きくなるにつれ、バーンイン（初期緩和）にかかるステップ数は大きくなる傾向がある。加えて初期緩和が起こった後の緩和もさらに遅くなり、古典ＭＣＭＣ計算では加速度的に状況が悪化する。対して量子試行ではこれが起こらず速やかにほぼ最適解に到達している。なお、両者共最適解に到達していないが、量子論で得られた解を古典ＭＣＭＣ計算の入力として用いれば、２都市の交換で最適解に到達できる。古典ＭＣＭＣベースの計算は、計算回数を１０倍に増やしても解の改善がほぼ起こらないため、計算が打ち切られている。古典ＭＣＭＣベースの計算では１０億回計算しても、量子試行の初期値までも下がっていない。

【0332】

＜＜量子―古典ハイブリッド計算＞＞
ulysses22を量子試行で解くことで得られた解である７０１９については、実質上問題が解けたとみなすことができる。これはスピン状態の微視的状態を示すことで分かる。ulysses22の場合、最適解は都市順序で並べると｛１，１４，１３，１２，７，６，１５，５，１１，９，１０，１９，２０，２１，１６，３，２，１７，２２，４，１８，８｝である。得られた解７０１３の微視的状態は都市１４と１３が入れ替わった状態であり、｛１，１３，１４，１２，７，６，１５，５，１１，９，１０，１９，２０，２１，１６，３，２，１７，２２，４，１８，８｝である。

【0333】

これは得られた解を他の解法の初期値として用いた場合、最低１回で真の解にたどり着けることを意味する。ＱＴＳＡによる解法はＭＣＭＣと異なるメカニズムで解候補を出す。そのためＱＴＳＡであれば、量子論だけで解にたどり着けない場合でも、古典論と組み合わせることで真の解に容易にたどり着ける場合がある。ulysses22の計算例は、古典論と組み合わせることで真の解に容易にたどり着ける一例である。

【0334】

ＱＴＳＡによる解法では、基底状態が確率的に求まるが、ＭＣＭＣの解法とは全く違い、量子力学に基づく一種のトンネル効果により解が求められる。古典ＭＣＭＣベースの解法では極小値に一旦はまると、それを抜け出すことは容易ではない。また、ハミング距離が大きく離れた状態にたどり着くのも容易ではない。しかし、量子論では関係が無い。熱的散逸からエネルギーが失われ、あたかもトンネル効果でポテンシャル障壁を透過するようにエネルギーが勝手に基底状態近傍に集まるからである。

【0335】

＜＜エネルギー分解能を定めることの効果＞＞
次に、式（５．３）から虚数時間発展法によるエネルギー分解能を予め決めてしまうことによる効果について説明する。エネルギー分解能は虚数時間発展のエネルギー収束判定の閾値を意味する。

【0336】

エネルギー分解能を定める効果を示すために実施した数値実験は次の通りである。問題は図１８に示した８都市のＴＳＰである。エネルギー分解能はΔＥ＝１０^-5、１０^-6、１０^-7、１０^-8、１０^-9、１０^-10、１０^-11、１０^-12とする。量子試行の回数はそれぞれ１０００回である。波動関数の初期状態に与えるノイズはランダムに与えられ、１粒子励起的にノイズが与えられる。ノイズの強度もランダムに決められている。このため、同じ量子ビットにノイズを与えた場合でも、異なる終状態に到達する。なぜならば、式（２．１６）に示す方程式は非線形であり、ノイズの振幅にも依存するからである。

【0337】

以下に数値計算の結果の概略を説明する。エネルギー分解能がΔＥ=１０^-6より大きい場合は古典に対応する状態に到達していない。これは不確定性原理の式（５．３）から、観測時間を短くした場合に対応する。量子アニーリングにおいて、高速に磁場を小さくするようなアニーリングをするとアニーリングに失敗することを意味する。

【0338】

ΔＥ=１０^-7の場合、１０００回の量子試行のうち８５８回が古典イジングモデルに対応する状態に収束し、残りの１４２回は古典論と対応する状態に収束していない。ΔＥ≦１０^-8の場合、１０００回全ての量子試行が全て古典論に対応する状態に収束している。

【0339】

これは量子アニーリングをする場合に磁場の消し方の速度に最大値があることを意味する。磁場の弱め方が速すぎると基底状態どころか古典論と対応する状態にも収束しない。
次にエネルギー分解能を小さくすることによる効果について説明する。

【0340】

図２９は、エネルギー分解能ごとの量子試行結果の一例を示す図である。図２９に示すグラフ７２は、エネルギー分解能をΔＥ＝１０^-7、およびΔＥ＝１０^-10にした場合の量子試行１０００回の結果得られたエネルギー値をヒストグラム表示したものである。グラフ７２の横軸はエネルギーであり、縦軸は頻度（該当エネルギーが得られた試行の回数）である。ヒストグラムを作るのに用いたエネルギー状態は古典論と対応する状態のみが用いられている。ΔＥ＝１０^-7の場合一定確率でサンプリングに失敗し、古典論と対応しない状態が得られるが、古典論と対応しないエネルギーの値はＥ≒１０²のオーダーである。

【0341】

このヒストグラムは量子試行において得られるエネルギーに関する確率分布を意味する。図２９からはエネルギー分解能を小さくすれば、エネルギー値が低い状態が多くサンプリングされることが分かる。また、基底状態をサンプリングする確率も上昇することが分かる。

【0342】

分布の形状はエネルギーの分解能が大きい時（粗い時）には、釣り鐘型になっているが、エネルギー分解能が小さくなるにつれ、基底状態のサンプリングに成功する確率があがるため、釣り鐘型の低エネルギー側の分布が０にならないような分布が形成されている。

【0343】

エネルギー分解能が小さくなるほど、低エネルギー側のサンプリングに成功する確率は上がるが、計算時間は伸びていく。虚数時間発展は指数関数的にエネルギーが減衰する傾向があるため、エネルギー分解能を０．１倍、０．０１倍としても計算時間は高々数倍の範囲に収まっている。しかし、エネルギー分解能を上げるほど計算時間は延びる。

【0344】

これはサンプリングにおいて最適な量子アニーリングのための横磁場の消し方の速度が存在することを意味する。横磁場を速く消せば計算時間は短くなるが、熱散逸機構が機能しづらくなるため、基底状態が得られる確率は下がる。逆に横磁場を十分遅くアニールすれば熱的散逸は機能するが、今度は計算時間が間延びする。

【0345】

なお、実験について言及しておくと、この計算の枠組みでは熱散逸効果を非ユニタリ項として取り込むが、途中で波動関数が０にならないように再規格化を行っている。もし実際の実験系に近づけたいのであれば、熱散逸項の部分に係数をかけ、熱散逸強度を調整することもできる。そして、その上で再規格化を行わないようにすれば、実際の系に近づくであろう。その場合、計算時間を長くすればするほど波動関数の振幅は消えてしまうことになり、最後に残るのはノイズだけになる。そのため、エネルギー分解能は実際の系では小さくしすぎることはできない。問題に応じて変わる最適な観測時間が存在する。つまり、アニーリング速度には最小値が存在する。

【0346】

波動関数の再規格化を行う枠組みでは、熱散逸により失われた確率が常に外部から増幅されるような理想的な状態を扱っている。
このように、１次の量子相転移が起きようと、系の基底状態を確率的に得られることを許容するなら、エネルギー分解能（アニーリング速度）は断熱遷移条件を守らなくてもよいことが分かる。

【0347】

〔第３の実施の形態〕
第２の実施の形態では、最適化装置１００は初期状態のスピンの波動関数にノイズを付加している。つまり、波動関数の対称性を破壊するような処置をすることで量子試行を定義している。しかし、ＱＴＳＡの支配方程式は平均場近似的である。ＱＴＳＡの解が確率的に基底状態に達するようにするためには、時間の進行と共に対称性が壊れればよい。そのため最適化装置１００は、磁場の係数に対してノイズを入れてもよい。

【0348】

例えば最適化装置１００は、ｉ番目のスピンの横磁場Ｊ_iに、式（８．１）のようにノイズを横磁場の係数に入れてもよい。

【0349】

【数75】

【0350】

ここで、Ｊは全体の横磁場の強度である。定式化の段階ではすべてのスピンが同じ横磁場強度を持っているとしたが、ノイズが加わるとして、最適化装置１００は横磁場の強度に揺らぎをランダムに入れる。このようにして、ほんの少しの大きさのノイズを付加するだけでも系の対称性を破壊することができる。そのため、波動関数にノイズを入れた場合と同様の効果が得られる。

【0351】

図３０は、横磁場にノイズを入れた場合のＱＴＳＡによる量子アニーリング結果の一例を示す図である。左側のグラフ８１は、横磁場の強さの関数としてエネルギーを示したものである。グラフ８１の横軸が横磁場の強さであり、縦軸がエネルギーである。グラフ８１において、黒丸印は横磁場にノイズを付加しない場合のＱＴＳＡの計算例である。白丸印と三角の印は共に、ノイズを付加した場合のＱＴＳＡの計算例である。

【0352】

右側のグラフ８２は、横磁場の関数として虚数時間発展が収束するまでに要した計算ステップ数を示したものである。グラフ８２の横軸が横磁場の強さであり、縦軸が虚数時間発展の収束に要した計算ステップ数である。黒丸印は横磁場にノイズを付加しない場合のＱＴＳＡの計算例である。白丸印はノイズを付加した場合のＱＴＳＡの計算例である。

【0353】

図３０に計算例として示したのは通常の意味での量子アニーリングプロトコルではエネルギーが下がらず、準安定状態から抜け出せない例である。グラフ８１に示すように、磁場に非対称性を入れて計算すると、横磁場が弱くなると、急激にエネルギーが下がるようになる。そして、横磁場を負の方向に大きくすると、今度は元の状態に戻ろうとする。つまり、何らかの形で対称性を壊すことにより系のエネルギーが下がる。

【0354】

またグラフ８２に示すように、非対称性を入れない場合（黒丸）にはほぼ一定の計算回数で虚数時間発展が収束する。Ｇ₂≒０の領域では計算回数が伸びるが、これは横磁場強度が弱くなるためスピンの反転が起こりづらくなるからである。半面、磁場に非対称性を入れると、複数のＧ₂の値で突発的に寿命が延びる。これは準安定状態が多数現れた状況である。このように非対称性をわずかに加えるだけで系の挙動は大きく変わり、よりエネルギーの低い状態へ崩壊していくことが分かる。

【0355】

このようにして磁場に対して非対称性を入れても、波動関数にノイズを入れた場合と同様の効果を得ることができる。なお、例示した磁場へのノイズの入れ方は一例である。波動関数に対するノイズ付加と同じように特定のｋ個のスピンを選んで一様なノイズを付加してもよいし、非一様なノイズを付加してもよい。

【0356】

〔その他の実施の形態〕
ノイズを付加したＱＴＳＡは、処理を並列化することも可能である。量子試行の定義から、ノイズを与えれば違う結果が得られる。そして、異なる２つの量子試行はお互いに独立に行える。従って異なるＬ台（Ｌは１以上の整数）の古典計算機に一度だけハミルトニアン情報を渡しておけば、あとは付加するノイズ情報だけ渡し、計算結果の量子ビット情報だけ回収すればよい。そして、演算もベクトル、行列演算だけである。このため、この計算はスーパースカラー型の並列計算機やＧＰＵなどのアクセラレータと相性が非常によい。

【0357】

また最適化装置１００は、イジング型で定義された最適化問題の係数を与えれば、ＱＴＳＡ法に基づき最適解を計算することができる。そのため最適化装置１００は、イジング型で定義できる最適化問題であれば、問題が属する分野や業種に依存せずに、最適解を計算可能である。例えば最適化装置１００は、金融、流通、材料工学、計算化学、計算生物学などの分野における最小値求解問題について求解可能である。また最適化装置１００は、金融、流通、ヘルスケア（製薬）などの業種における最小値求解問題について求解可能である。

【0358】

第２・第３の実施の形態では、プロセッサ１０１とメモリ１０２とを有するノイマン型コンピュータ構成によって最適化装置１００を実現しているが、例えば量子コンピュータを利用して最適化装置を実現することもできる。この場合、ノイマン型コンピュータによる最適化装置１００に、量子コンピュータによる最適化装置における熱的散逸を現象論的に取り込んでおいてもよい。これによりノイマン型コンピュータによる最適化装置１００を、量子コンピュータによる最適化装置のハードウェアの設計のシミュレータとしても利用可能となる。

【0359】

以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。

【符号の説明】

【0360】

１ イジングモデル
１０ 最適化装置
１１ 記憶部
１１ａ 求解問題情報
１２ 処理部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】

【図22】

【図23】

【図24】

【図25】

【図26】

【図27】

【図28】

【図29】

【図30】