特開 2025-102526 量子計算方法、量子計算システム、及びプログラム

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

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2025102526

(43)【公開日】2025-07-08

(54)【発明の名称】量子計算方法、量子計算システム、及びプログラム

(51)【国際特許分類】

   G06N 10/60 20220101AFI20250701BHJP

   G06N 99/00 20190101ALI20250701BHJP

【ＦＩ】

G06N10/60

G06N99/00 180

【審査請求】未請求

【請求項の数】5

【出願形態】ＯＬ

(21)【出願番号】P 2023220030

(22)【出願日】2023-12-26

【新規性喪失の例外の表示】特許法第３０条第２項適用申請有り Ａｄｉａｂａｔｉｃ Ｑｕａｎｔｕｍ Ｃｏｍｐｕｔｉｎｇ ２０２３、Ｔｈｅ Ｃｌｙｄｅ Ｈｏｔｅｌ（３３０ Ｔｉｊｅｒａｓ Ａｖｅ ＮＷ，Ａｌｂｕｑｕｅｒｑｕｅ ＮＭ ８７１０２，ＵＳ）、令和５年６月２１日（開催期間：令和５年６月１９日～令和５年６月２３日） ｈｔｔｐｓ：／／ａｒｘｉｖ．ｏｒｇ／ｓｕｂｍｉｔ／５３０９９７０、令和５年９月１５日

【国等の委託研究の成果に係る記載事項】（出願人による申告）令和元年度、国立研究開発法人科学技術振興機構、戦略的創造研究推進事業（ＣＲＥＳＴ）「地理空間情報を自在に操るイジング計算機の新展開」委託研究、産業技術力強化法第１７条の適用を受ける特許出願

(71)【出願人】

【識別番号】390001421

【氏名又は名称】学校法人早稲田大学

(74)【代理人】

【識別番号】110002675

【氏名又は名称】弁理士法人ドライト国際特許事務所

(72)【発明者】

【氏名】戸川 望

(72)【発明者】

【氏名】白井 達彦

(57)【要約】

【課題】コヒーレンス時間の限られた量子デバイスを用いて効率よく組合せ最適化問題を解く。
【解決手段】制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くために、量子デバイスは、イジングモデルの量子ゆらぎの強さを表すパラメータであるアニーリングパスにしたがって量子アニーリングを実行して、複数のスピン変数の量子状態を生成し（１０２）、量子状態の確率分布を生成する（１０４)。古典コンピュータは、量子状態を表す解に基づいて、制約を満たさない非実行可能解を制約を満たす実行可能解に変換する後処理を実行して（１０６）、後処理後の新たな確率分布を計算し（１０８）、新たな確率分布に基づいてイジングモデルのエネルギー期待値を計算し（１１０）、エネルギー期待値が小さくなるようにアニーリングパスを新たなアニーリングパスに更新する（１１２）。
【選択図】図１

【特許請求の範囲】

【請求項1】

制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算方法であって、
量子デバイスにより、前記イジングモデルの量子ゆらぎの強さを表すパラメータであるアニーリングパスにしたがって量子アニーリングを実行することにより、前記複数のスピン変数の量子状態を生成し、
前記量子デバイスにより、前記量子状態の確率分布を生成し、
古典コンピュータにより、前記量子状態を表す解に基づいて、前記制約を満たさない非実行可能解を前記制約を満たす実行可能解に変換する後処理を実行して、前記後処理後の新たな確率分布を計算し、
前記古典コンピュータにより、前記新たな確率分布に基づいて前記イジングモデルのエネルギー期待値を計算し、
前記古典コンピュータにより、前記エネルギー期待値が小さくなるように前記アニーリングパスを新たなアニーリングパスに更新する、量子計算方法。

【請求項2】

前記エネルギー期待値が収束するまで、前記量子状態の生成、前記確率分布の生成、前記後処理、前記エネルギー期待値の計算、及び前記アニーリングパスの更新のプロセスを繰り返すことにより、前記組合せ最適化問題の最適解を求める、請求項１に記載の量子計算方法。

【請求項3】

前記後処理は局所最適化手法に基づいて実行される、請求項１に記載の量子計算方法。

【請求項4】

制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算システムであって、
前記イジングモデルの量子ゆらぎの強さを表すパラメータであるアニーリングパスにしたがって量子アニーリングを実行することにより、前記複数のスピン変数の量子状態を生成し、前記量子状態の確率分布を生成する量子デバイスと、
前記量子状態を表す解に基づいて、前記制約を満たさない非実行可能解を前記制約を満たす実行可能解に変換する後処理を実行して、前記後処理後の新たな確率分布を計算し、前記新たな確率分布に基づいて前記イジングモデルのエネルギー期待値を計算し、前記エネルギー期待値が小さくなるように前記アニーリングパスを新たなアニーリングパスに更新する古典コンピュータと、
を備える、量子計算システム。

【請求項5】

制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算方法を量子デバイスと併用して古典コンピュータに実行させるためのプログラムであって、
前記量子デバイスによる量子アニーリングの実行により、前記複数のスピン変数の量子状態の確率分布が生成されると、前記量子状態を表す解に基づいて、前記制約を満たさない非実行可能解を前記制約を満たす実行可能解に変換する後処理を実行して、前記後処理後の新たな確率分布を計算し、
前記新たな確率分布に基づいて前記イジングモデルのエネルギー期待値を計算し、
前記エネルギー期待値が小さくなるように、前記量子アニーリングのパラメータであるアニーリングパスを新たなアニーリングパスに更新するステップを有するプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、制約のある組合せ最適化問題を解法するための量子計算方法、量子計算システム、及びプログラムに関する。

【背景技術】

【0002】

組合せ最適化問題とは、多数の組合せの中から最も良い組合せを選ぶために、制約条件のもとで目的関数を最小化又は最大化する変数の組合せを求める問題である。組合せ最適化問題は、一般に、量子アニーラやゲート型量子コンピュータなどの量子デバイス上でのイジングモデルの基底状態探索問題に変換される。近年、制約のある大規模な組合せ最適化問題を解法するための様々な量子アルゴリズムが提案されている（例えば、非特許文献１、非特許文献２、及び非特許文献３参照）。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】Edward Farhi, Jeffrey Goldstone, and Sam Gutmann, “A Quantum Approximate Optimization Algorithm,” arXiv:1411.4028 (2014).

【非特許文献2】Shunji Matsuura, Samantha Buck, Valentin Senicourt, and Arman Zaribafiyan, “Variationally scheduled quantum simulation,” Phys. Rev. A vol. 103, 052435 (2021).

【非特許文献3】Michiya Kuramata, Ryota Katsuki, and Kazuhide Nakata, “Larger sparse quadratic assignment problem optimization using quantum annealing and a bit-flip heuristic algorithm,” in 2021 IEEE 8th International Conference on Industrial Engineering and Applications (ICIEA), 2021, pp. 556-565.

【発明の概要】

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

【0004】

ノイズの影響により、量子デバイスがエラーなく量子アルゴリズムを実行することができる時間（コヒーレンス時間）には上限がある。しかしながら、従来の量子アルゴリズムを用いて大規模な組合せ最適化問題を解くためには非常に長い時間がかかるという課題がある。

【0005】

本発明は、上記課題に鑑みてなされたものであり、コヒーレンス時間の限られた量子デバイスを用いて効率よく組合せ最適化問題を解く量子計算方法、量子計算システム、及びプログラムを提供することを目的とする。

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

【0006】

本発明に係る量子計算方法は、制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算方法であって、量子デバイスにより、イジングモデルの量子ゆらぎの強さを表すパラメータであるアニーリングパスにしたがって量子アニーリングを実行することにより、複数のスピン変数の量子状態を生成し、量子デバイスにより、量子状態の確率分布を生成し、古典コンピュータにより、量子状態を表す解に基づいて、制約を満たさない非実行可能解を制約を満たす実行可能解に変換する後処理を実行して、後処理後の新たな確率分布を計算し、古典コンピュータにより、新たな確率分布に基づいてイジングモデルのエネルギー期待値を計算し、古典コンピュータにより、エネルギー期待値が小さくなるようにアニーリングパスを新たなアニーリングパスに更新する。

【0007】

本発明に係る量子計算システムは、制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算システムであって、イジングモデルの量子ゆらぎの強さを表すパラメータであるアニーリングパスにしたがって量子アニーリングを実行することにより、複数のスピン変数の量子状態を生成し、量子状態の確率分布を生成する量子デバイスと、量子状態を表す解に基づいて、制約を満たさない非実行可能解を制約を満たす実行可能解に変換する後処理を実行して、後処理後の新たな確率分布を計算し、新たな確率分布に基づいてイジングモデルのエネルギー期待値を計算し、エネルギー期待値が小さくなるようにアニーリングパスを新たなアニーリングパスに更新する古典コンピュータと、を備える。

【0008】

本発明に係るプログラムは、制約を有する組合せ最適化問題を複数のスピン変数で表されるイジングモデルを用いて解くための量子計算方法を量子デバイスと併用して古典コンピュータに実行させるためのプログラムであって、量子デバイスによる量子アニーリングの実行により、複数のスピン変数の量子状態の確率分布が生成されると、量子状態を表す解に基づいて、制約を満たさない非実行可能解を制約を満たす実行可能解に変換する後処理を実行して、後処理後の新たな確率分布を計算し、新たな確率分布に基づいてイジングモデルのエネルギー期待値を計算し、エネルギー期待値が小さくなるように、量子アニーリングのパラメータであるアニーリングパスを新たなアニーリングパスに更新するステップを有する。

【発明の効果】

【0009】

本発明によれば、量子デバイスがアニーリングパスにしたがって量子アニーリングを実行して量子状態の確率分布を生成した後、古典コンピュータが、非実行可能解を実行可能解に変換する後処理を実行し、後処理後の量子状態の新たな確率分布に基づいてエネルギー期待値が小さくなるようにアニーリングパスを更新するようにした。これにより、アニーリングパスが最適化され、コヒーレンス時間の限られた量子デバイスを用いて効率良く組合せ最適化問題を解くことが可能となる。

【図面の簡単な説明】

【0010】

【図1】本実施形態に係る量子計算方法を表すフローチャートである。

【図2A】２個のスピン変数で表される制約付き組合せ最適化問題において非実行可能解を実行可能解に変換する後処理を説明するための模式図である。

【図2B】後処理による確率分布の変化を説明するための模式図である。

【図3】３種類のアニーリングスケジュール（連続、線形、離散）を表す模式図である。

【図4】異なるアニーリング時間に対する連続スケジュールでの最適なアニーリングパスを表すグラフである。

【図5】複数の制約が互いに独立である例を説明するための模式図である。

【図6】複数の制約が互いに独立ではない例を説明するための模式図である。

【図7】複数の制約が互いに独立ではないときの後処理手法の一例を表す模式図である。

【図8】本実施形態に係る量子計算システムの構成を示すブロック図である。

【図9】図８に示す古典コンピュータのハードウェア構成を示すブロック図である。

【図10】グラフ分割問題の一例を表す模式図である。

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

【0011】

以下、図面を参照して、本発明の実施形態を説明する。

【0012】

まず、イジングモデルと量子アニーリングについて簡単に説明する。組合せ最適化問題（以下、ＣＯＰと呼ぶ。）は、一般に、イジングモデルの基底状態探索問題に変換される。イジングモデルは、無向グラフＧ＝（Ｖ、Ｅ）上で式（１）のように定義される。

【数1】

【0013】

ここで、Ｖは頂点の集合、Ｅは辺の集合である。σ_ｉ ^ｘ、σ_ｉ ^ｙ、σ_ｉ ^ｚは、それぞれ、サイトｉに作用するパウリ行列のｘ成分、ｙ成分、ｚ成分であり、Ｊ_ｉｊはスピン変数σ_ｉ ^ｚとσ_ｊ ^ｚとの間の相互作用係数を表し、ｈ_ｉはスピン変数σ_ｉ ^ｚ上の外部磁場係数を表し、Ｈ_０は定数である。

【0014】

制約付きＣＯＰのイジングモデルのハミルトニアンＨ_Isingは式（２）のように与えられる。Ｈ_obj、Ｈ_cstは、それぞれ、目的関数、制約項である。

【数2】

【0015】

制約を満たす解を実行可能解と呼び、制約を満たさない解を非実行可能解と呼ぶ。

【0016】

ここで、コスト関数Ｃ（σ）は式（３）のように与えられ、ｉ∈Ｖに対する｜σ_ｉ〉はσ_ｉ ^ｚの固有基底である。すなわち、σ_ｉ∈｛－１、１｝に対してσ_ｉ ^ｚ｜σ_ｉ〉＝σ_ｉ｜σ_ｉ〉である。

【数3】

【0017】

制約付きＣＯＰはquadratic unconstrained binary optimization（ＱＵＢＯ）モデルに変換して解くことができる。ＱＵＢＯモデルとは、イジングモデルを０と１の二値をとるスピン変数ｘ_ｉの二次式のエネルギー関数で表したものである。無向グラフＧ_ｑ＝（Ｖ_ｑ、Ｅ_ｑ）上で定義されたＱＵＢＯのハミルトニアンは式（４）のように与えられる。

【数4】

【0018】

ＱＵＢＯにおけるスピン変数ｘ_ｉ∈｛０、１｝を（σ_ｉ ^ｚ＋１）／２に置き換えると、式（１）で定義されるイジングモデルが得られる。

【0019】

量子アニーリングは、量子ゆらぎを用いてイジングモデルの基底状態を探索する手法である。量子アニーリングでは、式（５）で与えられるように、量子状態｜ψ（ｔ）〉が時間依存のシュレーディンガー方程式にしたがって時間発展する。ここで、プランク定数ｈをｈ／２π＝１としている。

【数5】

【0020】

式（５）におけるＨ（ｔ）は、イジングモデルのハミルトニアンＨ_Isingと量子ゆらぎに対応する横磁場ハミルトニアンＨ_ｑとにより、式（６）のように表される。

【数6】

【0021】

ｓ（ｔ）は時刻ｔでの量子ゆらぎの強さを表すパラメータであり、アニーリングパスと呼ばれる。ｓ（ｔ）は、量子アニーリングを実行する量子デバイスのハードウェアの制約により０と１との間の値をとる。アニーリングパスｓ（ｔ）のスケジューリングについては後述する（図３及び図４参照）。

【0022】

式（１）に示すように、イジングハミルトニアンＨ_Isingはパウリ行列のｚ成分σ_ｉ ^ｚで表されるが、横磁場ハミルトニアンＨ_ｑは、式（７）のように、パウリ行列のｘ成分σ_ｉ ^ｘで表される。

【数7】

【0023】

初期状態｜ψ（０）〉は横磁場ハミルトニアンＨ_ｑの基底状態に設定されている。アニーリング時間をＴと表記すると、出力状態｜ψ（Ｔ）〉は、ｔ＝０からｔ＝Ｔまで式（５）を解くことによって得られる。

【0024】

本実施形態に係る量子計算方法では、変分量子アルゴリズム（ＶＱＡ）と後処理とを用いて制約付きＣＯＰの最適解（基底状態）を探索する。ＶＱＡは、変分法によって短時間で（準）最適解に導く最適なアニーリングパスを見つける手法である。後処理は、量子デバイスの出力解をＣＯＰの制約を満たす実行可能解に変換する手法である。本実施形態では、ＶＱＡと後処理とを併用した手法を後処理変分量子アルゴリズム（ｐＶＱＡ）と呼ぶ。

【0025】

図１は、本実施形態の量子計算方法を表すフローチャートである。まず、ＣＯＰのイジングモデルとｔ∈［０、Ｔ］におけるアニーリングパスの初期値ｓ^(１)（ｔ）が量子デバイスに入力される。その後、量子デバイス及び古典コンピュータによってステップ１０２～１１２を繰り返すことでアニーリングパスを更新する。ｍ回目の量子アニーリングで用いるアニーリングパスをｓ^(m)（ｔ）、ｍ回の更新後に得られるアニーリングパスをｓ^(m+1)（ｔ）と表記する。

【0026】

量子デバイスは、現在のアニーリングパスｓ^(m)（ｔ）に沿って量子アニーリングを実行することによって量子状態｜ψ（Ｔ）〉を生成する（ステップ１０２）。

【0027】

次に、量子デバイスは、ステップ１０２で生成された量子状態｜ψ（Ｔ）〉の確率分布ｐ′（σ）を生成する（ステップ１０４）。ｐ′（σ）は式（８）のように定義される。実際には、量子状態｜ψ（Ｔ）〉の測定を繰り返すことで確率分布が近似的に求められる。

【数8】

【0028】

次に、古典コンピュータが、非実行可能解を実行可能解に変換する後処理を実行し（ステップ１０６）、後処理後の新たな確率分布ｐ（σ）を計算する（ステップ１０８）。後処理は、ステップ１０２で生成された量子状態｜ψ（Ｔ）〉を｛σ_ｉ ^ｚ｝基底で測定した際に得られた解σ′のマッピングＰ：Ｐσ′→σで定義される。後処理の具体例については後述する。新たな確率分布ｐ（σ）は式（９）のように与えられる。

【数9】

【0029】

次いで、古典コンピュータは、ステップ１０８で求めた新たな確率分布ｐ（σ）に基づいてエネルギー期待値Ｅを計算する（ステップ１１０）。エネルギー期待値Ｅは、コスト関数Ｃ（σ）（式（３））と確率分布ｐ（σ）（式（９））との加重和で表され、式（１０）のように定義される。

【数10】

【0030】

次いで、古典コンピュータは、ステップ１１０で計算されたエネルギー期待値Ｅが小さくなるように、現在のアニーリングパスｓ^(m)（ｔ）を新たなアニーリングパスｓ^(m+1)（ｔ）に更新する（ステップ１１２）。アニーリングパスの更新は、勾配降下法又はＰｏｗｅｌｌ法などの公知の最適化アルゴリズムを用いて行われる。

【0031】

ステップ１０２～１１２は、エネルギー期待値Ｅが収束するまで（解が更新されなくなるまで）繰り返される。上述のように、ステップ１０２及び１０４は量子デバイスによって実行され、ステップ１０６、１０８、１１０及び１１２は古典コンピュータによって実行される。

【0032】

ｐＶＱＡの適用例として、２個のスピン変数ｘ_１、ｘ_２∈｛０、１｝に対して、制約ｘ_１＋ｘ_２＝１の下で－ｘ_１を最小にするＣＯＰを考える。このＣＯＰのＱＵＢＯモデルにおけるハミルトニアン（エネルギー関数）は、例えば、式（１１）のように定式化される。ここで、Ｈ_obj＝－ｘ_１、Ｈ_cst＝５（ｘ_１＋ｘ_２－１）^２である。

【数11】

【0033】

図２Ａに、式（１１）で表されるＣＯＰの解のエネルギー値（Ｈの値）を示す。図２Ａにおいて、（ｘ_１、ｘ_２）＝（１、０）、（０、１）が実行可能解であり、（ｘ_１、ｘ_２）＝（０、０）、（１、１）が非実行可能解であり、（ｘ_１、ｘ_２）＝（１、０）が最適解である。

【0034】

後処理は、局所最適化手法に基づいて行うことができる。局所最適化手法に基づく後処理は、全ての極小解（local minimum solution）が実行可能解であるという性質を用いている。極小解は、スピン変数ｘ_ｉをフリップすることによってエネルギー値を下げない解として定義される。フリップとは、スピン変数の値を０から１又は１から０に変換することである。局所最適化手法は、スピン変数のフリップによってエネルギー値を下げる手法である。非実行可能解は極小解ではないので、局所最適化手法を繰り返してエネルギー値を下げていくと、極小解である実行可能解を得ることができる。

【0035】

局所最適化手法の例として貪欲法（greedy method）がある。図２Ａに示す制約付きＣＯＰに対し、貪欲法を用いて後処理する例を説明する。

【0036】

古典コンピュータは、暫定解（ｘ_１、ｘ_２）のエネルギー値Ｅ０と、（ｘ_１、ｘ_２）のｘ_１をフリップした解のエネルギー値Ｅ１と、（ｘ_１、ｘ_２）のｘ_２をフリップした解のエネルギー値Ｅ２とを計算する。Ｅ０、Ｅ１、Ｅ２の中でＥ１が最小であるとき、古典コンピュータは（ｘ_１、ｘ_２）のｘ_１をフリップする後処理を行う。Ｅ０、Ｅ１、Ｅ２の中でＥ２が最小であるとき、古典コンピュータは（ｘ_１、ｘ_２）のｘ_２をフリップする後処理を行う。Ｅ０、Ｅ１、Ｅ２の中でＥ０が最小であるときは、暫定解（ｘ_１、ｘ_２）が維持される。

【0037】

図２Ａの例では、貪欲法に基づく後処理によって、非実行可能解（０、０）、（１、１）はどちらも実行可能解（１、０）に変換される。一方、実行可能解である（１、０）、（０、１）は変換されない。よって、後処理によって実行可能解の確率が高くなる。

【0038】

図２Ｂに示すように、元の確率分布ｐ′は、非実行可能解（０、０）、（１、１）にも広がっている一方で、後処理後の新たな確率分布ｐは、実行可能解（０、０）、（１、１）のみに分布している。

【0039】

＜アニーリングスケジュール＞
次に、図３及び図４を参照して、ｐＶＱＡのパラメータであるアニーリングパスｓ（ｔ）のスケジューリングについて説明する。図３に示すように、アニーリングパスｓ（ｔ）には、連続スケジュール、線形スケジュール、離散スケジュールの３種類のスケジュールがある。

【0040】

連続スケジュールは、アニーリングパスｓ（ｔ）を連続関数で一般的に定式化したものであり、３種類のスケジュールの中で最も優れた性能を示すものの、既存の量子デバイスでの実現は困難を伴う。一方、線形スケジュールは量子アニーラで実行可能であり、離散スケジュールはゲート型量子コンピュータで容易に実行可能である。

【0041】

線形スケジュールは２つの変分パラメータｓ_１及びｓ_２を有し、それぞれ、ｔ＝０、ｔ＝Ｔでのｓ（ｔ）の値を表している。線形スケジュールのアニーリングパスｓ（ｔ）は式（１２）のように与えられる。

【数12】

【0042】

なお、従来の研究では、ｓ_１＝ｓ（０）＝０、ｓ_２＝ｓ（Ｔ）＝１に固定されていたが、本実施形態では、ｓ_１及びｓ_２を変分パラメータとして採用していることに注意されたい。

【0043】

離散スケジュールは量子近似最適化アルゴリズム（ＱＡＯＡ）に適用される。離散スケジュールは、ｐ個の層を有しており、各層は２個の変分パラメータで構成される。よって、２ｐ個の変分パラメータｓ_ｉ（ｉ∈｛１、…、２ｐ｝）が存在する。各層ｌ（ｌ∈｛１、…、ｐ｝）のアニーリングパスｓ（ｔ）は、式（１３）で表されるように０又は１の値をとる。

【数13】

【0044】

ここで、ｓ_０＝０であり、ｓ_ｉ－１≦ｓ_ｉ（ｉ∈｛１、…、２ｐ｝）、ｓ_２ｐ≦Ｔである。ｔ∈［ｓ_２ｐ、Ｔ］の間、量子状態は時間発展しない。式（６）及び式（１３）より、各層ｌにおいて量子状態は、イジングハミルトニアンＨ_Isingによる時間発展をした後に、横磁場ハミルトニアンＨ_ｑによる時間発展をする。

【0045】

図４に、異なるアニーリング時間（Ｔ＝０．１、Ｔ＝１、Ｔ＝１０）に対する連続スケジュールでの最適なアニーリングパスを示す。図４の各グラフにおける複数の線は、８個のノードを有する１０個の異なるグラフ分割問題に対するアニーリングパスに対応している。

【0046】

Ｔ＝０．１では、量子アニーリングはイジングハミルトニアンＨ_Isingで始まり、アニーリング時間Ｔの中間点で急に横磁場ハミルトニアンＨ_ｑに遷移する。Ｔ＝１では、量子アニーリングはイジングハミルトニアンＨ_Isingで始まり、連続関数が続き、横磁場ハミルトニアンＨ_ｑで終了する。Ｔ＝１０では、量子アニーリングは横磁場ハミルトニアンＨ_ｑで始まりイジングハミルトニアンＨ_Isingで終了しており、標準的なアニーリングパスに類似している。アニーリング時間が長くなると、上述のグラフ分割問題の最適解を得る成功確率は高くなる。ここで、成功確率は、図４に示す最適なアニーリングパスｓ（ｔ）に沿った量子アニーリングを計算し、その後、上述の後処理を実行したときに、最適解が得られる確率を示している。

【0047】

図２Ａでは、１個の制約を有するＣＯＰの例を示したが、本実施形態のｐＶＱＡは複数の制約を有するＣＯＰにも適用可能である。制約が１個のＣＯＰの例として、グラフ分割問題、ナップサック問題がある。互いに独立な複数の制約を有するＣＯＰの例として、グラフ彩色問題がある。互いに独立ではない複数の制約を有するＣＯＰの例として、巡回セールスマン問題、二次割当問題がある。

【0048】

図５に、複数の制約が互いに独立である例を示す。３行３列に並べられたスピン変数ｘ_ｉに対し、各行（横方向）のスピン変数ｘ_ｉの和が１という３個の制約があるものとする。この場合、制約は行ごとに独立に定義されるので、これら３個の制約は互いに独立である。図５の例では、１行目と２行目は制約を満たしているが、３行目はスピン変数ｘ_ｉの和が２となり、制約を満たしていない。

【0049】

図６に、複数の制約が互いに独立ではない例を示す。３行３列に並べられたスピン変数ｘ_ｉに対し、各行（横方向）のスピン変数ｘ_ｉの和が１、各列（縦方向）のスピン変数ｘ_ｉの和が１という６個の制約があるものとする。図６から明らかなように、１個のスピン変数ｘ_ｉに対して行と列の２個の制約が課せられていることから、これら６個の制約は互いに独立ではない。図６の例では、１行目、２行目、１列目、２列目は制約を満たしているが、３行目、３列目は制約を満たしていない。

【0050】

＜ｐＶＱＡの適用範囲＞
次に、ＣＯＰに関し、ｐＶＱＡの適用範囲について検討する。ＱＵＢＯモデルのハミルトニアンＱを式（１４－１）のように表す。ここで、Ｑ_objは目的関数、Ｑ_cstは制約項、Ａ＞０は制約係数である。さらに、後処理のためのハミルトニアンＱ′を式（１４－２）のように表す。ここで、Ｑ′_cst ^cは制約項、Ａ′＞０は制約係数である。

【数14】

【0051】

ｘ_ｉ∈｛０、１｝（ｉ∈Ｖ）に対し、式（１５）に示す線形の制約の下で目的関数Ｑ_objを最小化又は最大化するＣＯＰを考える。ここで、Ｍは制約の個数である。

【数15】

【0052】

このＣＯＰに関し、以下の定理を導くことができる。
定理：制約係数Ａ′の値が十分大きく（Ａ′＞ΔＱ_obj）、以下の第１条件、第２条件、及び第３条件を満たすとき、全ての非実行可能解を実行可能解に変換する局所最適化手法が存在する。ここで、ΔＱ_objは、スピン変数ｘ_ｉをフリップしたときのＱ_objのエネルギー変化の最大値である。
（ｉ）第１条件：式（１６）を満たすこと。

【数16】

（ｉｉ）第２条件：実行可能解が存在すること。
（ｉｉｉ）第３条件：式（１７）を満たすこと。すなわち、制約が互いに独立であること。

【数17】

【0053】

この定理を証明するためには、上述の条件の下で全ての非実行可能解が極小解ではないことを示す必要がある。ｘ_ｉ∈｛０、１｝（ｉ∈Ｖ）がｃ番目の制約を満足しないと仮定すると、式（１８）又は式（１９）が満たされる。

【数18】

【数19】

【0054】

式（１８）を満たすとき、第１条件から式（２０）が導かれる。

【数20】

【0055】

ｃ番目の制約を満足しない非実行可能解が極小解であると仮定して、背理法を用いた証明を行う。式（２０）と第３条件（すなわち、制約が互いに独立であること）は、スピン変数ｘ_ｉのフリップによって式（１５）のΣａ_ｉ ^cｘ_ｉの値が増加する場合、Ａ′＞ΔＱ_objのときエネルギー値（Ｑ′の値）が低くなることを示唆している。しかしながら、この非実行可能解が極小解であるならば、フリップによってエネルギー値は低くならないはずである。よって、フリップによってΣａ_ｉ ^cｘ_ｉの値は増加すべきではない。これは、ｘ_ｉ＝１のときａ_ｉ ^c≧０、ｘ_ｉ＝０のときａ_ｉ ^c≦０であることを意味し、任意のｘ_ｉ′に対して式（２１）が成り立つ。

【数21】

【0056】

しかしながら、式（２１）は、実行可能解が存在するという第２条件と矛盾している。よって、非実行可能解は極小解ではない。なお、式（１９）を満たす非実行可能解についても、同様の議論により、非実行可能解が極小解ではないことを証明することができる。

【0057】

複数の制約が互いに独立ではない場合、上述の定理に当てはまらず、局所最適化手法に基づく後処理を用いることができないが、アドホックな後処理方法を用いれば、ｐＶＱＡが適用可能である。図６に示したように、各行のスピン変数ｘ_ｉの和が１、各列のスピン変数ｘ_ｉの和が１という互いに独立ではない６個の制約を例に挙げてアドホックな後処理方法を説明する。

【0058】

スピン変数ｘ_ｉが、図７の初期状態７０２にあるとする。初期状態７０２は、２列目、３列目、１～３行目が制約を満たしていない。まず、ｘ_ｉの和が２以上の２列目及び３列目に対し、ｘ_ｉの和が１となり、且つ、エネルギー値ができる限り小さくなるようにｘ_ｉ＝１をｘ_ｉ＝０にフリップする。これにより、状態７０２から状態７０４に遷移したとする。

【0059】

状態７０４は、全ての列が制約を満たすが、２行目及び３行目が制約を満たしていない。次に、ｘ_ｉの和が２以上の２行目に対し、ｘ_ｉの和が１となり、且つ、エネルギー値ができる限り小さくなるようにｘ_ｉ＝１をｘ_ｉ＝０にフリップする。これにより、状態７０４から状態７０６に遷移したとする。

【0060】

状態７０６では、２列目、３列目、１行目、及び２行目が制約を満たしているが、１列目及び３行目ではｘ_ｉの和が０であり、制約を満たしていない。次に、ｘ_ｉの和が０の３行目に対し、ｘ_ｉの和が１となり、且つ、エネルギー値ができる限り小さくなるようにｘ_ｉ＝０をｘ_ｉ＝１にフリップする。これにより、状態７０６から状態７０８に遷移すれば、全ての制約が満たされる。このようなアドホックな後処理方法は、以下の非特許文献に記載されている。
S. Kanamaru, et al., “Mapping Constrained Slot-Placement Problems to Ising Models and its Evaluations by an Ising Machine,”2019 IEEE 9th International Conference on Consumer Electronics (ICCE-Berlin), Berlin, Germany, 2019.

【0061】

＜制約付きＣＯＰの例＞
次に、制約が１個のＣＯＰの例として、グラフ分割問題（ＧＰＰ）及び二次ナップサック（ＱＫＰ）を説明する。どちらも、ＮＰ困難なＣＯＰである。

【0062】

ＧＰＰは、偶数個の頂点（ノード）の集合Ｖ_ｇと辺集合Ｅ_ｇを有するグラフＧ_ｇ＝（Ｖ_ｇ、Ｅ_ｇ）が与えられ、Ｖ_ｇを等しいサイズの２つの部分集合に分割するとき、２つの部分集合を結ぶ辺の数（カット数）が最小になるように分割するものである。グラフの各頂点ｉ∈Ｖ_ｇに対してスピン変数ｘ_ｉを定義し、頂点ｉが第１部分集合に属するときｘ_ｉ＝１、頂点ｉが第２部分集合に属するときｘ_ｉ＝０と設定する。

【0063】

このＧＰＰは、式（２２）で与えられる線形の等式制約を有する。ここで、｜Ｖ_ｇ｜は頂点の数を表す。

【数22】

【0064】

このＧＰＰにおいて、式（１４－１）のＱ_obj及びＱ_cstは式（２３）のように与えられる。ここで、ｋ_ｉはノードｉの次数である。

【数23】

【0065】

ＱＫＰは、容量Ｃのナップサックとｎ個の品物が与えられたとき、容量Ｃを超えない範囲で、ナップサックに入れた品物の価値の合計を最大にする品物の組合せを見つけるＣＯＰである。ｉ番目の品物の容量及び価値を、それぞれ、ｗ_ｉ、ｐ_ｉｉと表記し、ｉ番目の品物とｊ番目の品物を同時にナップサックに入れたときの付加価値をｐ_ｉｊと表記する。ｉ番目の品物をナップサックに入れたときｘ_ｉ＝１と表し、入れないときｘ_ｉ＝０と表す。

【0066】

このＱＫＰは、式（２４）で与えられる線形の不等式制約を有する。

【数24】

【0067】

このＱＫＰにおいて、式（１４－１）のＱ_obj及びＱ_cstは式（２５）のように与えられる。

【数25】

【0068】

＜システム構成＞
次に、本実施形態の量子計算方法を実行する量子計算システムの構成について説明する。図８に示すように、本実施形態の量子計算システム１０は、古典コンピュータ２０と、古典コンピュータ２０に接続された量子デバイス３０とを備える。

【0069】

量子デバイス３０は、量子アニーラやゲート型量子コンピュータなどの量子デバイスである。量子デバイス３０は、図１のステップ１０２及び１０４に示したように、アニーリングパスに沿って量子アニーリングを実行して量子状態を生成し、量子状態の確率分布を生成する。

【0070】

古典コンピュータ２０は、パーソナルコンピュータなどのフォンノイマン型のコンピュータである。図９に示すように、古典コンピュータ２０は、プロセッサ２０２と、メモリ２０４と、記憶装置２０６と、入力部２０８と、ディスプレイ２１０と、I／Ｆ２１２とを備え、これらのデバイスはバスを介して接続される。

【0071】

プロセッサ２０２は、Central Processing Unit（ＣＰＵ）等を有し、メモリ２０４及び記憶装置２０６に格納されたプログラムにしたがって各種の処理を実行する。プロセッサ２０２は、図１のステップ１０６～１１２に示したように、後処理、新たな確率分布の計算、エネルギー期待値の計算、及びアニーリングパスの更新を行う。

【0072】

なお、プロセッサ２０２として、ＣＰＵ等の汎用コンピュータの代わりに、本実施形態の量子計算方法を実行するためのApplication Specific Integrated Circuits（ＡＳＩＣ）、Field Programmable Gate Array（ＦＰＧＡ）等の専用コンピュータを採用してもよい。

【0073】

メモリ２０４は、Read Only Memory（ＲＯＭ）及びRandom Access Memory（ＲＡＭ）を有する。ＲＯＭは、ＢＩＯＳ等のブートプログラムを格納している。プロセッサ２０２が、ＲＯＭに格納されたプログラム又は記憶装置２０６に格納されたプログラムを読み出す際、これらのプログラムはＲＡＭにロードされる。

【0074】

記憶装置２０６は、非一時的なコンピュータ読み取り可能な記憶媒体を有し、本実施形態の量子計算方法を量子デバイス３０と併用してプロセッサ２０２に実行させるためのプログラムと、プログラムの実行に必要なデータ等を格納している。コンピュータ読み取り可能な記憶媒体として、Hard Disk Drive（ＨＤＤ）、Solid State Drive（ＳＳＤ）、光ディスク等が挙げられる。

【0075】

なお、プロセッサ２０２によって実行されるプログラムを、ネットワークを介して接続された他のコンピュータに格納し、プロセッサ２０２が、他のコンピュータからI／Ｆ２１２を介してプログラムを読み出すようにしてもよい。

【0076】

入力部２０８は、マウス、キーボード等の入力デバイスを有する。ディスプレイ２１０は、Liquid Crystal Display（ＬＣＤ）等のディスプレイであり、プロセッサ２０２により実行された処理の結果を表示する。

【0077】

I／Ｆ２１２は、古典コンピュータ２０をLocal Area Network（ＬＡＮ）、Wide Area Network（ＷＡＮ）、及び／又はインターネット等のネットワークに接続するためのインターフェースである。

【0078】

次に、量子アニーラ及びゲート型量子コンピュータによってｐＶＱＡを実行した実施例１及び実施例２について説明する。実施例１ではＧＰＰを対象とし、実施例２ではＱＫＰを対象とする。

【実施例0079】

表１に、ＧＰＰの１０個の問題インスタンスＧＰＰ_ｎ_ｉに対してｐＶＱＡと従来の量子アニーリング（ＱＡ）とを量子アニーラで実行して得られた結果（残留エネルギーＲｅ及び成功確率ｐ_suc）を示す。ここで、ｎ＝｜Ｖ_ｇ｜、ｉはランダムグラフのラベルを表す。残留エネルギーＲｅは、コスト関数（式（３））の平均値と最適値との差を表す。成功確率ｐ_sucは、最適解が得られた確率を表す。コスト関数の平均値及び成功確率ｐ_sucは、量子アニーリングを繰り返すことで得られた１００個の解の中からコスト関数に基づく上位５０個の解に基づいて計算されたものである。線形のアニーリングパスを用い、グリッドサーチ法によりアニーリングパスを更新した。表１において、括弧内の数値は標準偏差を表している。Ｒｅの値が小さくｐ_sucの値が大きいほど、アルゴリズムの性能が良いことを示している。

【0080】

【表1】

【0081】

表１より、｜Ｖ_ｇ｜＝３２の５個のＧＰＰでは、ｐＶＱＡはＱＡに比べて、残留エネルギーＲｅが削減され、成功確率ｐ_sucが改善されていることがわかる。

【0082】

表２に、４ノードのＧＰＰ（｜Ｖ_ｇ｜＝４）に対してｐＶＱＡとＶＱＡとをゲート型量子コンピュータで実行して得られた成功確率ｐ_sucを示す。表２において、括弧内の数値は標準偏差を表している。図１０に、実際に用いられた４ノードのＧＰＰの模式図を示す。図１０における点線は最適解を表している。離散スケジュールのアニーリングパスを用い、Ｐｏｗｅｌｌ法によりアニーリングパスを更新した。アニーリングパスに沿った量子アニーリングを繰り返すことで得られた１００個の解から成功確率ｐ_sucを計算した。

【0083】

【表2】

【0084】

表２より、ｐＶＱＡはＶＱＡに比べて成功確率ｐ_sucが大きく、ｐＶＱＡがＶＱＡより著しく性能が優れていることがわかる。

【実施例0085】

表３に、ＱＫＰの７個の問題インスタンスＱＫＰ_ｎ_ｉに対してｐＶＱＡとＱＡとを量子アニーラで実行して得られた結果（残留エネルギーＲｅ及び成功確率ｐ_suc）を示す。ｎは品物の数、ｉはベンチマークであるインスタンスのラベルを表す。残留エネルギーＲｅ及び成功確率ｐ_sucの計算方法は表１と同様である。表３においても、括弧内の数値は標準偏差を表している。

【0086】

【表3】

【0087】

表３より、ｐＶＱＡはＱＡに比べて、残留エネルギーＲｅが大幅に削減され、成功確率ｐ_sucが改善されており、ｐＶＱＡがＱＡより著しく性能が優れていることがわかる。

【0088】

本実施形態の量子計算方法及びシステムによれば、アニーリングパスの変分最適化手法と後処理手法とを併用したアルゴリズムであるｐＶＱＡを用いることにより、コヒーレンス時間の限られた量子デバイスで効率良くＣＯＰを解くことができる。また、従来手法として解の精度を改善させることが可能となる。また、本実施形態の量子計算方法は、量子デバイスの種類を問わず適用可能である。

【0089】

なお、本発明は、上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲内で種々の変更が可能であり、当業者によってなされる他の実施形態、変形例も本発明に含まれる。