特開2020-204928(P2020-204928A)IP Force 特許公報掲載プロジェクト 2015.5.11 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 富士通株式会社の特許一覧
<>
  • 特開2020204928-最適化装置および最適化方法 図000029
  • 特開2020204928-最適化装置および最適化方法 図000030
  • 特開2020204928-最適化装置および最適化方法 図000031
  • 特開2020204928-最適化装置および最適化方法 図000032
  • 特開2020204928-最適化装置および最適化方法 図000033
  • 特開2020204928-最適化装置および最適化方法 図000034
  • 特開2020204928-最適化装置および最適化方法 図000035
  • 特開2020204928-最適化装置および最適化方法 図000036
  • 特開2020204928-最適化装置および最適化方法 図000037
  • 特開2020204928-最適化装置および最適化方法 図000038
  • 特開2020204928-最適化装置および最適化方法 図000039
  • 特開2020204928-最適化装置および最適化方法 図000040
  • 特開2020204928-最適化装置および最適化方法 図000041
  • 特開2020204928-最適化装置および最適化方法 図000042
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】特開2020-204928(P2020-204928A)
(43)【公開日】2020年12月24日
(54)【発明の名称】最適化装置および最適化方法
(51)【国際特許分類】
   G06N 99/00 20190101AFI20201127BHJP
【FI】
   G06N99/00 180
【審査請求】未請求
【請求項の数】10
【出願形態】OL
【全頁数】30
(21)【出願番号】特願2019-112547(P2019-112547)
(22)【出願日】2019年6月18日
(71)【出願人】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】110002918
【氏名又は名称】特許業務法人扶桑国際特許事務所
(72)【発明者】
【氏名】田村 泰孝
(72)【発明者】
【氏名】此島 真喜子
(57)【要約】
【課題】不等式制約を含む問題の求解を効率化すること。
【解決手段】エネルギー変化計算部12は、状態保持部11に保持された複数の状態変数の何れかの値が変化する場合に、複数の状態変数の値と重み値とに基づいて、複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算する。ペナルティ加算部13は、複数のエネルギーの変化値のそれぞれに、不等式制約における複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される、不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出する。更新制御部14は、設定された温度値と乱数値と複数の総エネルギーの変化値とに基づいて、状態保持部11に保持される複数の状態変数の何れかの値を変化させる。
【選択図】図1
【特許請求の範囲】
【請求項1】
エネルギーを表す評価関数に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する状態保持部と、
前記複数の状態変数の何れかの値が変化する場合に、前記複数の状態変数の値と前記重み値とに基づいて、前記複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算するエネルギー変化計算部と、
前記複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における前記複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される、前記不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出するペナルティ加算部と、
設定された温度値と乱数値と複数の前記総エネルギーの変化値とに基づいて、前記状態保持部に保持される前記複数の状態変数の何れかの値を変化させる更新制御部と、
を有する最適化装置。
【請求項2】
前記ペナルティ加算部は、前記接続係数と前記閾値とに基づいて算出される前記超過量の変化値に前記不等式制約の重要度を示す重要度係数を乗じることで前記ペナルティ値を算出する、
請求項1記載の最適化装置。
【請求項3】
前記不等式制約が上限制約値に対する不等式制約である場合、
前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第1の接続係数と第1の閾値とに基づいて算出される、前記第1の閾値に対する正方向の第1の超過量の変化値に応じた第1のペナルティ値を加算する、
請求項1または2記載の最適化装置。
【請求項4】
前記不等式制約が下限制約値に対する不等式制約である場合、
前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第2の接続係数と第2の閾値とに基づいて算出される、前記第2の閾値に対する負方向の第2の超過量の変化値に応じた第2のペナルティ値を加算する、
請求項1または2記載の最適化装置。
【請求項5】
前記不等式制約が上限制約値および下限制約値に対する不等式制約である場合、
前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第3の接続係数と第3の閾値とに基づいて算出される、前記第3の閾値に対する正方向の第3の超過量の変化値、および、前記第3の接続係数と前記第3の閾値よりも小さい第4の閾値とに基づいて算出される、前記第4の閾値に対する負方向の第4の超過量の変化値のうち、大きい方の変化値に応じた第3のペナルティ値を加算する、
請求項1または2記載の最適化装置。
【請求項6】
所定の制約値に対する等式制約の場合、
前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第4の接続係数と第5の閾値とに基づいて算出される、前記第5の閾値に対する正方向の第5の超過量の変化値、および、前記第5の閾値に対する負方向の第6の超過量の変化値のうち、大きい方の変化値に応じた第4のペナルティ値を加算する、
請求項1または2記載の最適化装置。
【請求項7】
前記状態保持部は、状態変数と前記不等式制約に対応する不等式制約用変数との組に対する前記接続係数を保持し、
前記接続係数は、前記不等式制約用変数の局所場における前記状態変数の重みを0以外の値とし、前記状態変数の局所場における前記不等式制約用変数の重みを0とする非対称接続係数である、
請求項1乃至6の何れか1項に記載の最適化装置。
【請求項8】
前記更新制御部により選択された変化対象の状態変数の識別情報に基づいて、前記エネルギー変化計算部により出力される前記複数のエネルギーの変化値のうちの何れかの変化値を選択し、選択した変化値を累積することで前記状態保持部に保持される前記複数の状態変数の値に対応するエネルギーの値を計算するエネルギー計算部を更に有する、
請求項1乃至7の何れか1項に記載の最適化装置。
【請求項9】
それぞれが前記状態保持部と前記エネルギー変化計算部と前記ペナルティ加算部と前記更新制御部とを備え、互いに異なる前記温度値が設定される複数の探索部と、
前記複数の探索部のそれぞれによる基底状態探索の繰り返し回数到達または一定時間経過後に、前記複数の探索部間で前記温度値または前記複数の状態変数の値を入れ替える交換制御部と、
を更に有する請求項1乃至8の何れか1項に記載の最適化装置。
【請求項10】
エネルギー変化計算部が、状態保持部に保持される複数の状態変数の何れかの値が変化する場合に、前記複数の状態変数の値と前記状態保持部に保持される重み値とに基づいて、前記複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算し、
ペナルティ加算部が、前記複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における前記複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される、前記不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出し、
更新制御部が、設定された温度値と乱数値と複数の前記総エネルギーの変化値とに基づいて、前記状態保持部に保持される前記複数の状態変数の何れかの値を変化させる、
最適化方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は最適化装置および最適化方法に関する。
【背景技術】
【0002】
ノイマン型コンピュータが不得意とする多変数の最適化問題を解く方法として、イジング型のエネルギー関数(評価関数や目的関数とも呼ばれる)を用いた最適化装置(イジングマシンやボルツマンマシンとも呼ばれる)がある。最適化装置は、計算対象の問題を、磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えて計算する。
【0003】
最適化装置は、例えば、ニューラルネットワークを用いてモデル化することもできる。その場合、イジングモデルに含まれる複数のスピンに対応した複数の状態変数のそれぞれが、他の状態変数の値と、他の状態変数と自身の状態変数との相互作用の大きさを示す重み係数とに応じて0または1を出力するニューロンとして機能する。最適化装置は、例えば、シミュレーテッド・アニーリングなどの確率的探索法により、上記のようなエネルギー関数の値(以下、エネルギーと言う)の最小値が得られる各状態変数の値の組合せを、解として求める。
【0004】
ここで、例えば、デジタル回路を用いてシミュレーテッド・アニーリングを行うことでエネルギーが最小となる各状態変数の値の組合せを計算する最適化装置の提案がある(例えば、特許文献1参照)。また、イジングモデルにおける各スピン間の相互作用を模擬するイジングチップとイジングチップを制御する情報処理部とを有する情報処理装置の提案もある(例えば、特許文献2参照)。
【0005】
なお、二値二次計画問題を半正定計画問題に緩和して、半正定計画問題の解を導出し、導出された半正定計画問題の解を二値二次計画問題の解に変換することで最適解を導出する最適化システムの提案もある(例えば、特許文献3参照)。
【先行技術文献】
【特許文献】
【0006】
【特許文献1】特開2018−41351号公報
【特許文献2】国際公開第2017/017807号
【特許文献3】国際公開第2017/056366号
【発明の概要】
【発明が解決しようとする課題】
【0007】
上記の最適化装置は、状態変数の二次形式で表される評価関数を最小化する。一方、最適化問題では、ある量(例えば、積載量や予算など)に対して上限や下限を設けるといった、不等式で表される制約条件(不等式制約)が課されることがある。しかし、こうした不等式制約を状態変数の二次形式で定式化した問題を、上記の最適化装置で解くのが難しいことがある。
【0008】
1つの側面では、本発明は、不等式制約を含む問題の求解を効率化する最適化装置および最適化方法を提供することを目的とする。
【課題を解決するための手段】
【0009】
1つの態様では、最適化装置が提供される。最適化装置は、状態保持部とエネルギー計算部とペナルティ加算部と更新制御部とを有する。状態保持部は、エネルギーを表す評価関数に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する。エネルギー変化計算部は、複数の状態変数の何れかの値が変化する場合に、複数の状態変数の値と重み値とに基づいて、複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算する。ペナルティ加算部は、複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される、不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出する。更新制御部は、設定された温度値と乱数値と複数の総エネルギーの変化値とに基づいて、状態保持部に保持される複数の状態変数の何れかの値を変化させる。
【0010】
また、1つの態様では、最適化方法が提供される。
【発明の効果】
【0011】
1つの側面では、不等式制約を含む問題の求解を効率化できる。
【図面の簡単な説明】
【0012】
図1】第1の実施の形態の最適化装置を示す図である。
図2】関数E(x)、C(x)の例を示す図である。
図3】変数間の関係の例を示す図である。
図4】最適化装置の回路構成例を示す図である。
図5】最適化装置の回路のブロック構成例を示す図である。
図6】判定部の回路構成例を示す図である。
図7】選択部の回路構成例を示す図である。
図8】変数毎のΔEおよびΔCの計算例を示す図である。
図9】関数E(x)に対応するエネルギーを計算する回路構成例を示す図である。
図10】最適化装置の演算例を示すフローチャートである。
図11】求解結果の比較の例を示す図である。
図12】最適化システムの例を示す図である。
図13】最適化装置の他の回路構成例を示す図である。
図14】第2の実施の形態の最適化装置の例を示す図である。
【発明を実施するための形態】
【0013】
以下、本実施の形態について図面を参照して説明する。
[第1の実施の形態]
第1の実施の形態を説明する。
【0014】
図1は、第1の実施の形態の最適化装置を示す図である。
最適化装置10は、計算対象の問題を変換したイジングモデルに含まれる複数のスピンに対応する複数の状態変数のそれぞれの値の組合せ(状態)のうち、エネルギー関数が最小値となるときの各状態変数の値(基底状態)を探索する。状態変数は、「バイナリ変数」や「ビット(スピンビット)」と呼ばれてもよい。
【0015】
最適化問題が不等式制約を含まない場合、最適化装置10は、次のように、エネルギー関数E(x)に基づいて基底状態の探索を行える。イジング型のエネルギー関数E(x)は、例えば以下の式(1)で定義される。
【0016】
【数1】
【0017】
式(1)に基づく制約なしの状態変数(バイナリ変数)の二次形式の関数の最適化問題は、(QUBO:Quadratic Unconstraint Binary Optimization)と呼ばれることがある。また、この二次形式での定式化は、QUBO形式と呼ばれることがある。
【0018】
式(1)の右辺第1項は、全状態変数から選択可能な2つの状態変数の全組合せについて、漏れと重複なく、2つの状態変数の値と重み値との積を積算したものである。xは、i番目の状態変数である。xは、j番目の状態変数である。重み値Wijは、i番目の状態変数とj番目の状態変数との間の重み(例えば、結合の強さ)を示す。なお、行列W={Wij}について、Wij=Wji、Wii=0であることが多い。状態変数xなどの変数に付加される添え字iは、当該変数の識別情報であり、インデックスと呼ばれる。
【0019】
式(1)の右辺第2項は、全状態変数のそれぞれのバイアス値と状態変数の値との積の総和である。bは、i番目の状態変数に対するバイアス値を示す。
例えば、イジングモデルにおけるスピンの「−1」は、状態変数の値「0」に対応する。イジングモデルにおけるスピンの「+1」は、状態変数の値「1」に対応する。
【0020】
状態変数xの値が変化して1−xとなると、状態変数xの増加分は、δx=(1−x)−x=1−2xと表せる。したがって、エネルギー関数E(x)に対して、状態変数xのスピン反転(値の変化)に伴うエネルギー変化ΔEは、式(2)で表される。
【0021】
【数2】
【0022】
は局所場(ローカルフィールド)と呼ばれ、式(3)で表される。
【0023】
【数3】
【0024】
状態変数xが変化(ビット反転)したときの局所場hの変化分δh(j)は、式(4)で表される。
【0025】
【数4】
【0026】
最適化装置10は、局所場hを保存するレジスタ(図示を省略している)を有し、状態変数xの値が変化したときに変化分δh(j)を、hに加算することで、ビット反転後の状態に対応するhを得る。
【0027】
最適化装置10では、基底状態の探索において、エネルギー変化がΔEとなる状態遷移(状態変数xの値の変化)を許容するか否かを決定するためにメトロポリス法やギブス法が用いられる。すなわち、最適化装置10は、ある状態から当該状態よりもエネルギーの低い他の状態への遷移を探索する近傍探索において、エネルギーが下がる状態だけでなく、エネルギーが上がる状態への遷移を確率的に許容する。例えば、エネルギー変化ΔEの状態変数の値の変化を受け入れる確率(遷移受入確率)Aは、式(5)で表される。
【0028】
【数5】
【0029】
ここで、βは温度Tの逆数(β=1/T)である。min演算子は、引数のうちの最小値を取ることを示す。例えば、メトロポリス法を用いる場合、式(5)におけるA=exp(−β・ΔE)に対して、両辺の自然対数を取って変形すると式(6)を得る。
【0030】
【数6】
【0031】
したがって、最適化装置10は、一様乱数u(0<u≦1)に対して、エネルギー変化ΔEが式(7)を満たす場合に、該当の状態変数の値の変化を許容する。
【0032】
【数7】
【0033】
そして、最適化装置10は、変化が許容された何れかの状態変数の値を変化させる。最適化装置10は、温度Tを初期温度から最低の温度まで下げながら、各温度において状態変数の値を変化させる処理を繰り返し実行することで、最適化問題に対する解を求める。
【0034】
最適化問題には不等式制約が付加されることがある。そこで、最適化装置10は、不等式制約に対する演算機能を提供する。最適化装置10は、状態保持部11、エネルギー変化計算部12、ペナルティ加算部13、更新制御部14、制約超過量計算部15および制約識別符号入力部16を有する。最適化装置10は、例えば、FPGA(Field Programmable Gate Array)などの半導体集積回路を用いて実現される。最適化装置10は、半導体チップとして実現され得る。
【0035】
状態保持部11は、エネルギーを表す評価関数(上記のエネルギー関数E(x))に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する。図1では、複数の状態変数の値を、状態ベクトルxで表す。例えば、複数の状態変数の数は、n個(nは2以上の整数)である(すなわち、状態ベクトルxはnビットである)。
【0036】
エネルギー変化計算部12は、状態保持部11に保持される複数の状態変数の何れかの値が変化されると、複数の状態変数の値と重み値とに基づいて、複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算する。例えば、エネルギー変化計算部12は、変化候補を状態変数xとする場合のエネルギーの変化値ΔEを、j=1〜nのそれぞれに対して求める。
【0037】
ペナルティ加算部13は、複数の状態変数に対して計算された複数のエネルギーの変化値ΔEのそれぞれに、不等式制約に違反する超過量(制約超過量)に応じたペナルティ値ΔCを加算して総エネルギーの変化値(ΔE+ΔC)を算出する。制約超過量は、不等式制約における複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される値である。例えば、変化候補の状態変数xに対応する制約超過量の変化値が、制約超過量計算部15により計算され、ペナルティ加算部13に供給される。ペナルティ加算部13は、状態変数xに対応する制約超過量の変化値に基づいてペナルティ値ΔCを求め、ΔEに加算する。例えば、ペナルティ加算部13は、当該制約超過量の変化値に、不等式制約に対する重要度を表す所定の係数を乗じることでペナルティ値ΔCを算出することが考えられる。
【0038】
更新制御部14は、設定された温度値と乱数値と複数の総エネルギーの変化値とに基づいて、状態保持部11に保持される複数の状態変数の何れかの値を変化させる。更新制御部14は、複数の総エネルギーの変化値と熱励起エネルギー(熱ノイズ)との相対関係によって複数の状態変数の何れかを受け入れるか否かを確率的に決定するともいえる。更新制御部14は、遷移可否判定部14aおよび選択部14bを有する。
【0039】
遷移可否判定部14aは、温度値Tと乱数値uと総エネルギーの変化値(ΔE+ΔC)とに基づいて、当該総エネルギーの変化値(ΔE+ΔC)に対応する状態遷移を許容するか否かを示すフラグFを出力する。例えば、状態変数xの値の変化による状態遷移を許容する場合にF=1であり、状態変数xの値の変化による状態遷移を許容しない場合にF=0である。なお、温度値Tは、図示を省略している制御部(あるいは、制御回路)により、遷移可否判定部14aに入力される。ここで、遷移可否判定部14aによる判定は、式(7)のΔEを、総エネルギーの変化値(ΔE+ΔC)に置き換えた式(8)に基づいて行われる。
【0040】
【数8】
【0041】
選択部14bは、遷移可否判定部14aにより出力されたF(j=1〜n)のうち、F=1に対応する状態遷移の番号(インデックス)を1つ選択する。選択部14bは、選択したインデックスを状態保持部11に出力することで、状態保持部11に保持される複数の状態変数の何れかの値を変化させる。例えば、状態保持部11は、選択部14bから受け付けたインデックスに対応する状態変数の値を変化(ビット反転)させる。こうして、状態保持部11に入力されるクロック信号(clk)に同期して、状態保持部11に保持される状態変数の値が更新される。
【0042】
制約超過量計算部15は、不等式制約における閾値と、不等式制約に対する状態変数xの接続係数とに基づいて、状態変数xを変化候補とした場合の制約超過量の変化値を計算し、ペナルティ加算部13に供給する。閾値は、不等式制約における上限値でもよいし、下限値でもよい。閾値は、上限値および下限値の両方でもよく、その場合、不等式制約に違反する制約超過量の変化値として、上限値に対する制約超過量の変化値および下限値に対する制約超過量の変化値のうちの大きい方を採用することが考えられる。
【0043】
制約識別符号入力部16は、制約識別符号を制約超過量計算部15に入力する。制約識別符号は、閾値として用いられる上限値や下限値を指定する情報である。制約識別符号は、前述の制御部により、制約識別符号入力部16に入力される。
【0044】
また、最適化装置10で保持可能な状態変数の最大数をN(Nはnより大きい整数)として、N個の変数のうちn個を通常の状態変数として用い、K(Kは1以上の整数)個を不等式制約を表すための変数(不等式制約用変数)として用いる構成が考えられる。N=n+Kである。Kは、不等式制約の数に対応し、不等式制約の数が1個であればK=1、不等式制約の数が2個であればK=2などとなる。
【0045】
その場合、制約識別符号入力部16は、N個の変数のうち、状態ベクトルに属する通常の状態変数と、状態ベクトルに属さない不等式制約用変数とを、制約識別符号により指定する。そのため、制約識別符号入力部16は、エネルギー変化計算部12、ペナルティ加算部13および遷移可否判定部14aに対して制約識別符号を入力する構成としてもよい(図1では当該構成を示している)。このようにすると、後述されるように、K個の不等式制約用変数のそれぞれに対応する局所場hを記憶するレジスタに格納された現在のhの値を用いてペナルティ値ΔCを計算することができる。
【0046】
N個の変数のうち、n個を通常の状態変数として用い、K個を不等式制約用変数として用いる場合、不等式制約が付加された最適化問題を次のように定式化できる。ここで、通常の状態変数のインデックスの集合をΨと表す。集合Ψの要素数はn個である。不等式制約用変数のインデックスの集合をΩ(Ω={k,k,…,k})と表す。集合Ωの要素数はK個である。
【0047】
最小化するエネルギー関数は、n個の状態変数の関数E(x)として次の式(9)で表される。
【0048】
【数9】
【0049】
n個の状態変数から構成されるK個の不等式制約は式(10)で表される。
【0050】
【数10】
【0051】
接続係数aijは、インデックスiに対応する不等式制約(不等式制約iと表記する)における状態変数x(j∈Ψ)の重みを示す。Cliは、不等式制約iにおける下限値である。Cuiは、不等式制約iにおける上限値である。ただし、不等式制約として、下限値または上限値の何れか一方のみを設けてもよい。
【0052】
このとき、総エネルギーに対応する目的関数Etot(x)を式(11)で表す。
【0053】
【数11】
【0054】
関数C(x)は、不等式制約に対する制約超過量と重要度係数λの積を、全ての不等式制約について和をとったものであり、式(12)で表される。
【0055】
【数12】
【0056】
ここで、重要度係数λは、0以上の実数であり、不等式制約iの重要度を表す。また、max演算子は、引数のうちの最大値をとることを示す。関数C(x)の変化値がペナルティ値に相当する。
【0057】
図2は、関数E(x)、C(x)の例を示す図である。
図2(A)は、関数E(x)の一例のグラフG1を示す。図2(B)は、関数C(x)の一例のグラフG2を示す。グラフG1,G2の横軸は、状態x(状態ベクトルで表される状態)を示す。グラフG1の縦軸はE(x)である。グラフG2の縦軸はC(x)である。関数E(x)に関数C(x)を加えた目的関数Etot(x)に関して、例えば、式(8)(メトロポリス法)に基づいて状態を遷移させることで、関数E(x)に対して不等式制約を充足する(あるいは不等式制約に関する制約超過を最小化する)解を得られる。
【0058】
図3は、変数間の関係の例を示す図である。
変数群g1は、状態変数の集合{x}(j∈Ψ)を示す。変数群g2は、不等式制約用変数の集合{x}(k∈Ω)を示す。
【0059】
状態変数のペアは、双方向のカップリングである。すなわち、当該ペアの何れか一方の状態変数の反転に対して、当該ペアの他方の状態変数の局所場が影響を受ける関係となる。
【0060】
状態変数と不等式制約用変数とのペアは、状態変数から不等式制約用変数へ向かう単方向のカップリングである。すなわち、状態変数の反転に対して、不等式制約用変数の局所場が影響を受けるが、不等式制約用変数の反転に対して、状態変数の局所場は影響を受けない関係となる。不等式制約kに対する状態変数xの影響の大きさ(重み)は、接続係数akjにより表される。接続係数akjは、前述の行列Wの要素として与えられる。また、不等式制約用変数のバイアスbは、b=0が設定される。
【0061】
状態変数と不等式制約用変数との関係を実装するため、一例では、最適化装置10は、不等式制約用変数xの値を0に設定し、不等式制約用変数の反転が生じないように制御することで、状態変数の局所場への影響を抑制する。他の例では、接続係数akjを非対称接続係数(すなわち、akj≠0,ajk=0)とすることで、不等式制約用変数xによる状態変数の局所場への影響を抑制し得る。更に他の例では、遷移可否判定部14aにより、制約識別符号に基づいて、不等式制約用変数xに対して常にフラグF=0を出力するように制御することも考えられる。
【0062】
不等式制約用変数は他の不等式制約用変数とカップリングされない。すなわち、不等式制約用変数のペアに対しては、Wij,Wji(i,j∈Ω)を0とする。
不等式制約用変数x(k∈Ω)の局所場hは式(13)で表される。
【0063】
【数13】
【0064】
不等式制約用変数は、通常の状態変数で表される状態のエネルギーE(x)には影響を与えないが、不等式制約用変数xに対応する局所場hの値は、通常の状態変数の反転に伴って、最新の値に更新される。そして、局所場hと、不等式制約の上限・下限とを比較することで、不等式制約の充足状況が分かる。
【0065】
式(11)で表されるように、総エネルギーEtot(x)は、QUBO項E(x)と不等式制約項(ペナルティ項)C(x)との和となる。通常の状態変数x(j∈Ψ)の反転に対する総エネルギーのエネルギー変化ΔEtot j(x)は、式(14)で表される。
【0066】
【数14】
【0067】
ΔEは、例えば、特開2018−41351号公報に開示された既存の回路構成を用いて、QUBO項に基づき生成される。
ΔCは、不等式制約用変数x(i∈Ω)の局所場hの現在の値を用いて、式(15)に基づき生成される。
【0068】
【数15】
【0069】
ここで、ΔCijは、式(16)により表される。
【0070】
【数16】
【0071】
ただし、δxは、状態変数x反転時の状態変数xの変化分であり、式(17)で表される。
【0072】
【数17】
【0073】
すなわち、ΔCijは、状態変数xの反転に伴うインデックスiの不等式制約(不等式制約i)に対する制約超過量の変化分の値(変化値)である。また、ペナルティ値ΔCは、不等式制約毎の制約超過量の変化値の、重要度係数による重み付き総和である。なお、不等式制約iで上限値のみを与える場合、式(15)、(16)のmax演算において、「Cli−aijδx−h」の引数や「Cli−h」の引数を含まなくてよい。また、不等式制約iで下限値のみを与える場合、式(15)、(16)のmax演算において、「h+aijδx−Cui」の引数や「h−Cui」の引数を含まなくてよい。更に、制約iとして等式制約を表す場合、式(15)、(16)において、Cui=Cliとすればよい。
【0074】
次に、最適化装置10に対する上記原理の実装例を更に説明する。以下では、最適化装置10で保持可能な状態変数の最大数をN個とする。また、Ωに属するインデックスについて、便宜的にkを1、kを2,…,kをKのように、kの添え字iで表記することがある。
【0075】
図4は、最適化装置の回路構成例を示す図である。
エネルギー変化計算部12は、ΔE計算部12a1〜12aNを有する。ΔE計算部12a1〜12aNのそれぞれは、状態保持部11から供給されるj番目の状態変数を反転候補とする場合のエネルギーの変化値ΔEを、式(2)〜(4)に基づいて計算し、ペナルティ加算部13に出力する。ΔE計算部12a1〜12aNのそれぞれには、局所場h〜hを記憶するレジスタ(図示を省略している)から、局所場h〜hの現在の値が供給される。
【0076】
なお、ΔE計算部12a1〜12aNには、不等式制約用変数に対応するΔE計算部が含まれるが、後述されるように、不等式制約用変数に対応するΔE計算部は、制約識別符号に基づいて、比較的大きな値を出力することで、不等式制約用変数の反転を抑制する。
【0077】
ペナルティ加算部13は、係数レジスタ13a1〜13aK、乗算器13b1〜13bK、加算器13c11〜13c1N,…,13cK1〜13cKNを有する。ここで、係数レジスタ13a1〜13aK、乗算器13b1〜13bKのように、不等式制約用変数のインデックスに対応するK個の係数レジスタおよび乗算器を図示しているが、当該係数レジスタおよび乗算器は、それぞれN個設けられる。すなわち、不等式制約用変数のインデックスに対応していない他の係数レジスタおよび他の乗算器の図示を省略している。同様に、加算器13c11〜13c1N,…,13cK1〜13cKNのように、不等式制約用変数に対応するK群(あるいはK行)の加算器を図示しているが、当該加算器群は、N群(あるいはN行)設けられる。すなわち、不等式制約用変数のインデックスに対応していない他の加算器群の図示を省略している。
【0078】
係数レジスタ13a1〜13aKは、式(12)、(15)における重要度係数λ〜λを保持し、係数レジスタ13a1〜13aKに対応する乗算器13b1〜13bKに供給する。
【0079】
乗算器13b1〜13bKは、制約超過量計算部15により計算されたΔC1j〜ΔCKj(j∈Ψ)に、係数λ〜λを乗じて、加算器13c11〜13c1N,…,13cK1〜13cKNに供給する。
【0080】
具体的には、乗算器13b1は、ΔC11にλを乗じて、加算器13c11に供給する。また、乗算器13b1は、ΔC12にλを乗じて、加算器13c12に供給する。以降、同様に、乗算器13b1は、ΔC1Nにλを乗じて、加算器13c1Nに供給する。
【0081】
更に、乗算器13bKは、ΔCK1にλを乗じて、加算器13cK1に供給する。また、乗算器13bKは、ΔCK2にλを乗じて、加算器13cK2に供給する。以降、同様に、乗算器13bKは、ΔCKNにλを乗じて、加算器13cKNに供給する。他の乗算器についても同様である。
【0082】
加算器13c11〜13c1N,…,13cK1〜13cKNは、乗算器13b1〜13bKから供給されるλΔC11〜λΔC1N,…,λΔCK1〜λΔCKNを、エネルギー変化計算部12で計算されたΔE〜ΔEに加算する。
【0083】
具体的には、加算器13c11は、ΔE計算部12a1から供給されるΔEにλΔC11を加算し、インデックス1に対応する次段の加算器(図示を省略している)に供給する。加算器13c12は、ΔE計算部12a2から供給されるΔEにλΔC12を加算し、インデックス2に対応する次段の加算器に供給する。以降、同様に、加算器13c1Nは、ΔE計算部12aNから供給されるΔEにλΔC1Nを加算し、インデックスNに対応する次段の加算器に供給する。
【0084】
更に、加算器13cK1は、インデックス1に対応する前段の加算器(図示を省略している)から供給される値にλΔCK1を加算し、遷移可否判定部14aに出力する。加算器13cK2は、インデックス2に対応する前段の加算器から供給される値にλΔCK2を加算し、遷移可否判定部14aに出力する。以降、同様に、加算器13cKNは、インデックスNに対応する前段の加算器から供給される値にλΔCKNを加算し、遷移可否判定部14aに出力する。
【0085】
遷移可否判定部14aは、判定部14a1〜14aNを有する。判定部14a1〜14aNのそれぞれは、ペナルティ加算部13から供給されるΔE+ΔC〜ΔE+ΔCに基づいて、状態変数xの反転を許容するか否かを、式(8)により確率的に判定し、反転可否を示すフラグF〜Fを選択部14bに供給する。判定部14a1〜14aNの回路構成の詳細は後述される。
【0086】
選択部14bは、遷移可否判定部14aから供給されるフラグF〜Fに基づいて、反転させる状態変数を選択し、当該状態変数のインデックスを状態保持部11に供給する。選択部14bの回路構成の詳細は後述される。
【0087】
制約超過量計算部15は、ΔC計算部15a1〜15aKを有する。ΔC計算部15a1〜15aKは、式(16)に基づいて、ΔC1j〜ΔCKj(j∈Ψ)を計算し、ペナルティ加算部13に供給する。
【0088】
具体的には、ΔC計算部15a1は、ΔC11〜ΔC1Nを、乗算器13b1に供給する。以降、同様に、ΔC計算部15aKは、ΔCK1〜ΔCKNを乗算器13bKに供給する。
【0089】
上記のように、最適化装置10の機能は、電子回路によって実現できる。例えば、状態保持部11、エネルギー変化計算部12、ペナルティ加算部13、遷移可否判定部14a、選択部14bおよび制約超過量計算部15を、それぞれ状態保持回路、エネルギー変化計算回路、ペナルティ加算回路、遷移可否判定回路、選択回路および制約超過量計算回路のように呼ぶこともできる。
【0090】
図5は、最適化装置の回路のブロック構成例を示す図である。
図4で例示した回路要素は、インデックスi=1,2,…,Nのそれぞれに対応するブロックBiに分けて実装することができる。図5では、最適化装置10のj番目のブロックBjの回路構成例が示されている。
【0091】
ブロックBjは、ΔE計算部12aj、係数レジスタ13a1j〜13aKj、乗算器13b1j〜13bKj、加算器13c1j〜13cKj、判定部14ajおよびΔC計算部15a1j〜15aKjを有する。
【0092】
ΔE計算部12ajは、インデックスjに対応する状態変数xの値、および、局所場hの値に基づいて、ΔEを計算し、加算器13c1jに供給する。
係数レジスタ13a1j〜13aKjは、係数λ〜λを保持する。係数レジスタ13a1j〜13aKjをブロック毎に設けるのではなく、各ブロックで共有のレジスタとしてもよい。
【0093】
乗算器13b1j〜13bKjは、それぞれ係数レジスタ13a1j〜13aKjに保持されるλ〜λを、ΔC計算部15a1j〜15aKjから供給されるΔC1j〜ΔCKjに乗算して、加算器13c1j〜13cKjに供給する。乗算器13b1j〜13bKjは、それぞれ乗算器13b1〜13bKのうちのインデックスjに対応する回路要素であると言える(例えば、乗算器13b1jは乗算器13b1のインデックスjに対応する回路である)。
【0094】
加算器13c1jは、ΔE計算部12ajから供給されるΔEに、乗算器13b1jから供給されるλΔC1jを加算して、次段の加算器(図示を省略している)に供給する。以降、同様にλkj(k∈Ω)が順次加算される。加算器13cKjは、前段の加算器(図示を省略している)から供給される値にλΔCKjを加算し、その結果(ΔE+ΔCに相当)を判定部14ajに供給する。
【0095】
判定部14ajは、総エネルギーの変化値ΔE+ΔCに基づいて、インデックスjに対応する反転可否を判定し、判定結果に応じたフラグFを選択部14bに出力する。
ΔC計算部15a1j〜15aKjは、それぞれΔC1j〜ΔCKjを計算し、乗算器13b1j〜13bKjに供給する。ΔC計算部15a1j〜15aKjは、それぞれΔC計算部15a1〜15aKのうちのインデックスjに対応する回路要素であると言える(例えば、ΔC計算部15a1jはΔC計算部15a1のインデックスjに対応する回路要素である)。
【0096】
ここで、ブロックBjには、インデックスjが通常の状態変数に対応するか、または、不等式制約用変数に対応するかを示す制約識別符号が、制約識別符号入力部16により入力される。制約識別符号は、制約無し(不等式制約用変数に該当しない)、上限のみの制約、下限のみの制約、あるいは、上限および下限の両方の制約かを識別する識別コードを含む。
【0097】
例えば、識別コードdは、2ビット(di1,di2)で表される。識別コードd=00は、インデックスiが通常の状態変数のインデックスであることを示す。識別コードd=01は、インデックスiが上限のみの制約(あるいは上限のみの制約に対応する不等式制約用変数)のインデックスであることを示す。識別コードd=10は、インデックスiが下限のみの制約(あるいは下限のみの制約に対応する不等式制約用変数)のインデックスであることを示す。識別コードd=11は、インデックスiが上限および下限の両方の制約(あるいは上限および下限の両方の制約に対応する不等式制約用変数)のインデックスであることを示す。
【0098】
制約識別符号は、識別コードに対する上限値および下限値の指定も含む。例えば、制約識別符号は、(識別コード,上限値,下限値)の組となる。
等式制約の場合は、識別コードd=11を用い、前述のように、上限値および下限値を一致させる。ΔE計算部12ajは、識別コードに応じて、インデックスjに対するΔEを計算する。ΔC計算部15a1j〜15aKjは、識別コードに応じて、インデックスjに対するΔC1j〜ΔCKjを計算する。識別コードに応じたΔEおよびΔC1j〜ΔCKjの計算方法の例は後述される。
【0099】
図6は、判定部の回路構成例を示す図である。
判定部14ajは、オフセット値生成部21、乱数生成部22、ノイズ値生成部23、符号反転回路24、加算器25,26および比較器27を有する。
【0100】
オフセット値生成部21は、選択部14bから出力される遷移可否を示すフラグに基づいて、オフセット値Eoff(Eoff≧0)を生成し、加算器25に供給する。具体的には、オフセット値生成部21は、選択部14bから出力されるフラグが遷移可を示す場合、オフセット値を0にリセットする。オフセット値生成部21は、選択部14bから出力されるフラグが遷移不可を示す場合、オフセット値に増分値ΔEoffを加算する。当該フラグが連続して遷移不可を示す場合、オフセット値生成部21は、ΔEoffを積算することで、EoffをΔEoffずつ増加させる。
【0101】
乱数生成部22は、0<u≦1の一様乱数uを生成し、ノイズ値生成部23に出力する。
ノイズ値生成部23は、使用法則(例えば、メトロポリス法)に応じた所定の変換テーブルを保持する。ノイズ値生成部23は、一様乱数uと、制御部(あるいは制御回路)により供給された温度Tを示す温度情報とを用いて、当該変換テーブルにより、式(8)に基づくノイズ値(熱ノイズ)に対応する−T・ln(u)の値を生成する。ノイズ値生成部23は、生成した−T・ln(u)の値を加算器26に出力する。
【0102】
符号反転回路24は、加算器13cKjから供給される総エネルギーの変化値(ΔE+ΔC)の符号を反転させ、加算器25に供給する。
加算器25は、符号反転回路24から供給される−(ΔE+ΔC)にオフセット値Eoffを加算し、加算器26に供給する。
【0103】
加算器26は、加算器25から供給される−(ΔE+ΔC)+Eoffに、熱ノイズ−T・ln(u)を加算し、比較器27に供給する。
比較器27は、加算器26により出力された評価値(−(ΔE+ΔC)+Eoff−T・ln(u))を閾値(具体的には0)と比較することで、式(8)に基づく判定を行う。比較器27は、評価値が0以上の場合、遷移可を示すフラグ(F=1)を、選択部14bに出力する。比較器27は、評価値が0未満の場合、遷移不可を示すフラグ(F=0)を、選択部14bに出力する。
【0104】
ここで、選択部14bから出力されるフラグが遷移不可を示す場合、現在の状態が局所解に陥っていると考えられる。オフセット値生成部21による、−(ΔE+ΔC)へのEoffの加算やEoffの漸増により、状態遷移が許容されやすくなり、現在の状態が局所解にある場合、その局所解からの脱出が促進される。
【0105】
次に、選択部14bの回路構成例を説明する。
図7は、選択部の回路構成例を示す図である。
選択部14bは、複数段にわたってツリー状に接続された複数のセレクタ部および乱数ビット生成部32a,32b,…,32rを有する。乱数ビット生成部32a〜32rは、ツリー状に接続された複数のセレクタ部の段毎に設けられる。乱数ビット生成部32a〜32rの各々は、0又は1の値をとる1ビット乱数を生成し、各段の選択回路に供給する。1ビット乱数は、入力されたフラグのペアのうちの何れか一方を選択するために用いられる。
【0106】
初段のセレクタ部31a1,31a2,…,31apの各々には、判定部14a1〜14aNの各々が出力する遷移可否を示すフラグの組が入力される。例えば、セレクタ部31a1には、1番目の判定部14a1が出力するフラグと、2番目の判定部14a2が出力するフラグとのペアが入力される。また、セレクタ部31a2には、3番目の判定部が出力するフラグと、4番目の判定部が出力するフラグとのペアが入力される。以降、同様にして、セレクタ部31apには、N−1番目の判定部が出力するフラグと、N番目の判定部14aNが出力するフラグとのペアが入力される。このように、隣り合う判定部が出力するフラグのペアが初段のセレクタ部に入力される。初段のセレクタ部31a1〜31apの数は、N/2となる。以降、段を経る毎にセレクタ部の数は半分になる。
【0107】
セレクタ部31a1〜31apの各々は、入力されたフラグのペアと、乱数ビット生成部32aが出力する1ビット乱数に基づいて、入力されたフラグのペアのうちの一方を選択する。セレクタ部31a1〜31apの各々は、選択したフラグと選択したフラグに対応する1ビットの識別値とを含む状態信号を、2段目のセレクタ部31b1〜31bqに出力する。例えば、セレクタ部31a1が出力する状態信号と、セレクタ部31a2が出力する状態信号とがセレクタ部31b1に入力される。同様に、セレクタ部31a1〜31apの隣り合うセレクタ部が出力する状態信号のペアが、2段目のセレクタ部に入力される。
【0108】
セレクタ部31b1〜31bqの各々は、入力された状態信号のペアと、乱数ビット生成部32bが出力する1ビット乱数に基づいて、入力された状態信号のうちの一方を選択する。セレクタ部31b1〜31bqの各々は、選択した状態信号を、3段目のセレクタ部31c1,…に出力する。ここで、セレクタ部31b1〜31bqは、選択した方の状態信号に含まれる状態信号について、何れの状態信号を選択したかを示すように1ビットを付加して更新し、選択した状態信号を出力する。
【0109】
3段目以降のセレクタ部においても同様の処理が行われ、各段のセレクタ部で1ビットずつ識別値のビット幅が増えていき、最終段のセレクタ部31rから、選択部14bの出力である状態信号が出力される。選択部14bが出力する状態信号に含まれる識別値が、2進数で表されたインデックスに対応する。ただし、図7で示す回路構成例では、変数のインデックスを0始まりとする場合としている(1始まりのインデックスに対応させるには、選択部14bが出力する識別値に1を加算した値をインデックスとすればよい)。
【0110】
例えば、図7では、セレクタ部31bqの回路構成例が示されている。2段目以降の他のセレクタ部も、セレクタ部31bqと同様の回路構成により実現される。
セレクタ部31bqの入力は、1つ目の状態信号(status_1)と、2つ目の状態信号(status_2)である。セレクタ部31b1の出力は、状態信号(status)である。セレクタ部31bqは、OR回路41、NAND回路42およびセレクタ43,44を有する。
【0111】
OR回路41には、状態信号(status_1)に含まれるフラグ(flag1)と、状態信号(status_2)に含まれるフラグ(flag2)とが入力される。例えば、状態信号(status_1)が前段の2つのセレクタ部のうちの上位側(indexの大きい側)の出力、状態信号(status_2)が前段の2つのセレクタ部のうちの下位側(indexの小さい側)の出力である。OR回路41は、flag1とflag2とのOR演算結果(flag)を出力する。
【0112】
NAND回路42には、flag1とflag2とが入力される。NAND回路42は、flag1とflag2とのNAND演算結果を、セレクタ43の選択信号入力端子に出力する。
【0113】
セレクタ43には、flag1と1ビット乱数(rand)とが入力される。セレクタ43は、NAND回路42から入力されるNAND演算結果に基づいて、flag1又はrandの何れかを選択し、出力する。例えば、セレクタ43は、NAND回路42のNAND演算結果が「1」の場合、flag1を選択し、NAND回路42のNAND演算結果が「0」の場合、randを選択する。
【0114】
セレクタ44には、状態信号(status_1)に含まれる識別値(index1)と、状態信号(status_2)に含まれる識別値(index2)とが入力される。セレクタ44の選択信号入力端子には、セレクタ43の選択結果が入力される。セレクタ44は、セレクタ43の選択結果に基づいて、index1又はindex2の何れかを選択し、出力する。例えば、セレクタ44は、セレクタ43の選択結果が「1」の場合、index1を選択し、セレクタ43の選択結果が「0」の場合、index2を選択する。
【0115】
OR回路41およびセレクタ43,44の出力の組が、セレクタ部31bqが出力する状態信号(status)となる。
なお、初段のセレクタ部の入力は識別値を含まない。このため、初段のセレクタ部は、識別値(図中ではindexと表記)として、選択した方に対応するビット値(下位側の場合に「0」、上位側の場合に「1」)を追加して出力する回路となる。
【0116】
このように、選択部14bは、遷移可であるスピンビットから1つをトーナメント方式で選択する。トーナメントの各試合(すなわち、各選択回路での選択)では、勝った方(すなわち、選択された方)のエントリ番号(0又は1)をindex wordの上位ビットに付け加えていく。最終段のセレクタ部31rが出力するインデックスが選ばれたスピンビットを示す。例えば、変数の数Nが1024個の場合、最終段のセレクタ部31rが出力する状態信号は、遷移可否を示すフラグと、10ビットで表されるインデックスとを含む。
【0117】
ただし、インデックスの出力方法は、上記のように選択部14bで生成する方法以外の方法も考えられる。例えば、判定部14a1〜14aNの各々から選択部14bに各判定部に対応するインデックスを供給し、選択部14bにより遷移可否を示すフラグとともに、当該フラグに対応するインデックスを選択してもよい。この場合、判定部14a1〜14aNのそれぞれは、自身に対応するインデックスを格納するインデックスレジスタを更に有し、インデックスレジスタから選択部14bにindexを供給する。
【0118】
図8は、変数毎のΔEおよびΔCの計算例を示す図である。
最適化装置10は、全てのi(i=1〜N)に関する(x,h),(x,h),…,(x,h)、および、識別コードdに対して下記の手順を並列に実行する。
【0119】
(S1)最適化装置10は、i∈Ωであるか否かを判定する。i∈Ωでない場合、ステップS2,S2aに処理が進む。i∈Ωである場合、ステップS3,S3aに処理が進む。i∈Ωであるか否かは、d=(di1,di2)により判定できる。すなわち、di1およびdi2の少なくとも一方が1(d=01,10,11)ならば、i∈Ωである。また、di1およびdi2の両方が0(d=00)ならば、i∈Ωでない。
【0120】
(S2)最適化装置10は、ΔEを式(2)により計算される値(ΔE=(2x−1)h)とする。
(S2a)最適化装置10は、全てのインデックスjにわたって、ΔCij=0とする。
【0121】
(S3)最適化装置10は、全てのインデックスjにわたって、行列Wを保持するレジスタR1に記憶された係数aijに基づき、ΔCijを式(16)により計算される値とする。ΔCijの計算に用いられる局所場hは、hを保持するレジスタから、不等式制約iに対するΔCijを計算する各ΔC計算部に供給される。
【0122】
(S3a)最適化装置10は、ΔE=Emaxとする。ここで、Emaxは、非常に大きな値であり、判定部14aiにおける判定で常にF=0となる値である。例えば、EmaxをΔEとしてとれる値の最大値とすることが考えられる。
【0123】
また、最適化装置10は、総エネルギーEtot(x)のうち、QUBO項に相当するエネルギーE(x)を次の回路によってモニタリングする。
図9は、関数E(x)に対応するエネルギーを計算する回路構成例を示す図である。
【0124】
最適化装置10は、更に、エネルギー計算部17および制御部18を有する。
エネルギー計算部17は、総エネルギーEtot(x)のうち、QUBO項に相当するエネルギーE(x)を計算する。エネルギー計算部17は、更新制御部14により選択された変化対象の状態変数の識別情報に基づいて、エネルギー変化計算部12により出力される複数のエネルギーの変化値のうちの何れかの変化値を選択する。エネルギー計算部17は、選択した変化値を累積することで状態保持部11に保持される複数の状態変数の値に対応するエネルギーの値E(x)を計算する。具体的には、エネルギー計算部17は、選択部17aおよび累算器17bを有する。
【0125】
選択部17aは、反転させる状態変数のインデックスiの選択結果を選択部14bから受け付ける。選択部17aは、エネルギー変化計算部12により出力されるΔE(i=1,…,n)のうち、選択部14bから供給されるインデックスiに対応するΔEを選択し、累算器17bに出力する。
【0126】
累算器17bは、エネルギーE(x)の値を保持し、選択部17aから供給されるΔEを、エネルギーE(x)の値に累算することで、現在の状態変数の値に対応するエネルギーE(x)を更新し、当該エネルギーE(x)を出力する。なお、エネルギーE(x)の初期値は、状態ベクトルの初期値に対して、演算開始時に累算器17bに与えられる。
【0127】
制御部18は、最適化装置10における探索部D1による演算を制御する。ここで、探索部D1は、状態保持部11、エネルギー変化計算部12、ペナルティ加算部13、更新制御部14、制約超過量計算部15、制約識別符号入力部16、選択部17aおよび累算器17bを含む回路ブロックである。制御部18は、探索部D1における探索回数(状態変数の反転回数)や温度設定などを制御する。制御部18は、電子回路により実現される。制御部18は、制御回路とも呼べる。
【0128】
最適化装置10は、図9の回路構成により、QUBO項に相当するエネルギーE(x)をモニタリングできる。図9の回路構成によってエネルギーE(x)をモニタリングすることで、不等式制約における不等式制約項における係数λの値が比較的大きい場合でも、エネルギーE(x)の計算の桁あふれがないという利点がある。
【0129】
次に、最適化装置10による演算の手順を説明する。
図10は、最適化装置の演算例を示すフローチャートである。
(S10)制御部18は、探索部D1の初期化を行う。例えば、制御部18は、温度や、温度毎の繰り返し回数C1および温度の更新回数C2などの外部からの設定を受け付け、探索部D1に初期温度を設定する。また、制御部18は、制約識別符号を制約識別符号入力部16に設定する。制御部18は、初期状態を状態保持部11に設定する。このとき、制御部18は、制約識別符号により不等式制約用変数として指定された変数の値を0に設定する。更に、制御部18は、重みW、係数λ、各変数に対応する局所場の初期値を所定のレジスタに設定するとともに、初期状態に対応するエネルギーE(x)を累算器17bに設定する。制御部18は、初期化が完了すると探索部D1による演算を開始させる。
【0130】
下記のステップS11,S12についてブロックBjに着目して説明するが、ステップS11,S12は、ブロックB1〜BNにおいて、並列に実行される。
(S11)ブロックBjは、hとxとを読み出す。前述のようにhは、所定のレジスタに保持される。xは、状態保持部11に保持される。
【0131】
(S12)ブロックBjは、図8の手順に基づいてΔCijおよびΔEの計算を行うとともに、ΔCijに基づいてΔCを計算し、ΔEtot j=ΔE+ΔCを求める。ブロックBjは、−(ΔE+ΔC)に対してEoffを加算し、更に、熱ノイズ−T・ln(u)を加算した評価値と閾値との判定に基づいて、フラグFを出力する。
【0132】
(S13)選択部14bは、ブロックB1〜BNにより出力されるフラグF〜Fに基づいて、反転対象のインデックスiを選択し、選択したインデックスiを状態保持部11に供給する。
【0133】
(S14)状態保持部11は、選択部14bから供給されるインデックスiに対応する状態変数xの値を反転させる(変化させる)ことで、状態を更新するとともに、所定のレジスタに保持された各変数に対応する局所場hを更新する(式(4)に基づく更新)。なお、本ステップの説明では、変化対象のインデックスをiと表記している。当該更新により、式(13)に基づく不等式制約用変数に対する局所場の値も更新される。
【0134】
(S15)制御部18は、現温度において、規定回数C1だけ状態変数の更新を行ったか否かを判定する。規定回数C1だけ更新を行った場合、ステップS16に処理が進む。規定回数C1だけ更新を行っていない場合、ステップS11に処理が進む。
【0135】
(S16)制御部18は、温度の更新を行う。これにより、探索部D1に設定される温度Tが下げられる。制御部18は、規定回数C1の更新のたびに、予め定められた方法により、探索部D1に設定する温度Tを下げる。
【0136】
(S17)制御部18は、規定回数C2だけ温度を更新したか否かを判定する。規定回数C2だけ温度を更新した場合、ステップS18に処理が進む。規定回数C2だけ温度を更新していない場合、ステップS11に処理が進む。
【0137】
(S18)制御部18は、累算器17bにより出力された最低エネルギー(エネルギーE(x)の最小値)に対応する状態を出力する。例えば、制御部18は、探索部D1が最終的に到達した状態を、最低エネルギーに対応する状態として出力してもよい。あるいは、制御部18は、上記の探索の過程で、探索部D1が到達した最低エネルギーに対応する状態を出力してもよい。
【0138】
ここで、ステップS12において、ペナルティ加算部13は、次のようにペナルティ値をエネルギーの変化値に加算することで、総エネルギーの変化値ΔEtot j=ΔE+ΔC(j∈Ψ)を算出する。
【0139】
例えば、ペナルティ加算部13は、接続係数aij(i∈Ω)と閾値(CuiおよびCli)に基づいて算出される制約超過量の変化値ΔCijに不等式制約iの重要度を示す重要度係数λを乗じることでペナルティ値ΔCを算出し、ΔEに加算する。ペナルティ加算部13は、不等式制約が複数ある場合は、不等式制約毎のλΔCijをiについて和をとって、ΔCを求める。重要度係数により、QUBO項に対する不等式制約の重みや、不等式制約が複数ある場合の各不等式制約の重みを調整可能にして、求解を行える。
【0140】
また、不等式制約が上限制約値Cui(i∈Ω)に対する不等式制約である場合、ペナルティ加算部13は、複数のエネルギーの変化値ΔEのそれぞれに対して、第1の接続係数aijと第1の閾値Cuiとに基づいて算出される、第1の閾値Cuiに対する正方向の第1の超過量の変化値に応じた第1のペナルティ値ΔCを加算する。
【0141】
また、不等式制約が下限制約値Cliに対する不等式制約である場合、ペナルティ加算部13は、複数のエネルギーの変化値ΔEのそれぞれに対して、第2の接続係数aijと第2の閾値Cliとに基づいて算出される、第2の閾値Cliに対する負方向の第2の超過量の変化値に応じた第2のペナルティ値ΔCを加算する。
【0142】
また、不等式制約が上限制約値Cuiおよび下限制約値Cliに対する不等式制約である場合、ペナルティ加算部13は、複数のエネルギーの変化値ΔEのそれぞれに対して、第3の接続係数aijと第3の閾値Cuiとに基づいて算出される、第3の閾値Cuiに対する正方向の第3の超過量の変化値、および、第3の接続係数aijと第3の閾値Cuiよりも小さい第4の閾値Cliとに基づいて算出される、第4の閾値Cliに対する負方向の第4の超過量の変化値のうち、大きい方の変化値に応じた第3のペナルティ値ΔCを加算する。
【0143】
また、所定の制約値に対する等式制約の場合、ペナルティ加算部13は、複数のエネルギーの変化値ΔEのそれぞれに対して、第4の接続係数aijと第5の閾値Cui=Cliとに基づいて算出される、第5の閾値Cui=Cliに対する正方向の第5の超過量の変化値、および、第5の閾値Cui=Cliに対する負方向の第6の超過量の変化値のうち、大きい方の変化値に応じた第4のペナルティ値ΔCを加算する。
【0144】
図11は、求解結果の比較の例を示す図である。
図11(A)は、不等式制約を含む最適化問題に対する最適化装置10を用いた求解結果の例を示す。図11(B)は、ペナルティ加算部13を有していない既存の最適化装置を用いた、不等式制約を状態変数の二次形式で表す場合(二乗制約の場合)の最適化問題に対する求解結果の例を示す。
【0145】
ここで、図11では、最適化問題の一例として、ナップザック問題を考える。(ペナルティ加算部13を有していない)既存の最適化装置を用いてナップザック問題を解く場合、次のように問題を定式化し得る。
【0146】
ここで考えるナップザック問題は、N個の品物(インデックスiの品物の重さはw、値段はv)を、積載容量がCのナップザックに入れ、ナップザック中の品物の値段の合計(総価値)を最大にする問題である。ナップザックに品物iを入れる場合に状態変数x=1とし、入れない場合に状態変数x=0とする。
【0147】
この場合、ナップザック中の品物の総価値Vは、式(18)で表される。
【0148】
【数18】
【0149】
また、ナップザックにおける容量制約は、式(19)で表される。
【0150】
【数19】
【0151】
式(19)の容量制約は、スラック変数S(S≧0)を用いて、式(20)のように表される。
【0152】
【数20】
【0153】
式(20)に基づいて、QUBOの目的関数Eを式(21)で表せる。
【0154】
【数21】
【0155】
Sを式(22)のようにバイナリ展開し、状態変数x,yのQUBO形式として解くことが考えられる。状態変数x,yのQUBO形式は式(23)で表される。
【0156】
【数22】
【0157】
【数23】
【0158】
スラック変数Sは、最適化で比較的小さな値となるため、S=0としても最適解が取得され得る。よって、式(24)を得る。
【0159】
【数24】
【0160】
そこで、図11の比較では、次のようなテストデータを作成して実験を行った。
ナップザック容量C=4952991である。品物毎の重さW=D〜Dは、それぞれ区間[0,100]の一様乱数に100000を加算した値である。品物毎の値段P=P〜Pは、それぞれ区間[1,1000]の一様乱数として発生させた値である。
【0161】
また、実験条件を次の通りとした。第1に、不等式制約に対する重要度係数λをλ=10^zとして、z=−10から10まで、1刻みで網羅する。第2に、温度T=5,10,50,100,500(初期温度)として、それぞれメトロポリス基準に従ったモンテカルロ法による演算を行う。第3に、乱数uの初期値を100個用意し、この中で1つでも制約を満たす解が得られる初期値があれば、「解が得られる」ものとする。
【0162】
図11の比較における最適化装置10に対する不等式制約(絶対値制約)を含むエネルギー関数E(x)は、式(25)で表される。
【0163】
【数25】
【0164】
また、図11の比較における二乗制約の場合のエネルギー関数E(x)は、式(26)で表される。
【0165】
【数26】
【0166】
ただし、式(26)では、スラック変数は小さな値が最適であるため、スラック変数を0と置いた。
図11(A)および図11(B)の各グラフの横軸はλ=10^zにおけるzである。図11(A)および図11(B)の各グラフの縦軸は価値Pの合計値(最大)である。なお、制約を満たさない場合には、P=0である。
【0167】
上記の実験条件では、絶対値制約では何れの温度でもλが10^(−1)以上であれば、解を得られるのに対して、二乗制約では全てP=0となり、制約を満たす解を得られなかった。そこで、λの範囲を広げ、温度TをT=10^10,10^15,10^20というように大きく設定したところ、二乗制約では限られた範囲のみで解が得られ、二乗制約では探索に比較的大きなコストがかかることが分かる。
【0168】
係数λが比較的小さい場合に解が得られない理由は、総価値に比べて不等式制約項の値が小さくなることで、制約が軽視され、制約を満たす解を得られる可能性が低下するためと考えられる。
【0169】
図11(A)において、系列p1は、T=10^10の場合である。系列p2は、T=10^15の場合である。系列p3は、T=10^20の場合である。系列p4は、T=5の場合である。系列p5は、T=10の場合である。系列p6は、T=50の場合である。系列p7は、T=100の場合である。系列p8は、T=500の場合である。
【0170】
図11(B)において、系列q1は、T=10^10の場合である。系列q2は、T=10^15の場合である。系列q3は、T=10^20の場合である。図11(B)の二乗制約の例では、温度T=5,10,50,100,500に対応する系列は、解が得られなかったため、図示を省略している。
【0171】
二乗制約の場合において、このような結果となった原因は、二乗制約のエネルギー関数の式にあると考えられる。すなわち、ナップザックに入っている品物の重さの合計がナップザックの容量付近である場合、品物を取り除いても式(26)に基づくエネルギーが増加してしまい、探索ができなくなるためである。特に、二乗制約の場合は、二次の項により、一次の項で表される絶対値制約の場合よりも、式(26)のエネルギー関数E(x)の値に対する不等式制約項の寄与が過大になり得る。
【0172】
ここで、最適化問題では様々な制約条件が入るため、二次形式だけでは問題を定式化するのか困難なことがある。上記の例以外にも、例えば、配車配送問題などの積載量上限、工場を建設するための予算の上限など、必要なリソース量が上限を超えない、あるいは下限を下回らないなど不等式で表される制約が頻出し得る。
【0173】
これらの不等式制約は、二次形式で表して既存の最適化装置で解くことが難しいことがあり、また、不等式制約を厳密に満たすように探索すると最適解に到達し難いという問題がある。これは、不等式制約を破るような中間状態を経由しないと最適解に到達できないことが多いためである。
【0174】
そこで、ペナルティ加算部13を有する最適化装置10を用いることで、不等式制約を二次形式に変形しなくても、元の一次式のままで容易に扱えるようになる。
図11で例示されるように、最適化装置10によれば、既存の最適化装置よりも、係数λや温度Tについて、広い範囲で解を得ることができる。このため、二乗制約の場合よりも、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。すなわち、最適化装置10の応用範囲を広げることができ、より多くの問題の求解に、最適化装置10を利用可能になる。
【0175】
上記の例では、k∈Ωに対して、x=0かつΔE=Emaxとすることで、変数xによるエネルギーE(x)への影響を抑制するものとしたが、係数akjとして、非対称接続係数を用いて、変数xによるエネルギーE(x)への影響を抑制してもよい。非対称接続係数は、前述のように、akj≠0、ajk=0と表される。なお、非対称接続係数は、変数間の重み値の記号Wを用いてWkj≠0、Wjk=0と表されてもよい。
【0176】
すなわち、状態保持部11は、状態変数と不等式制約に対応する不等式制約用変数との組に対する接続係数aijを保持してもよい。接続係数aijは、不等式制約用変数x(k∈Ω)の局所場hにおける状態変数xの重み(akj)を0以外の値とし、状態変数xの局所場hにおける不等式制約用変数xの重み(ajk)を0とする非対称接続係数であってもよい。
【0177】
非対称接続係数を用いる場合、選択部14bにより不等式制約用変数が反転対象として選択されたとしても、エネルギーE(x)に対する変数xの反転の影響を抑制できる。
次に、第1の実施の形態の他の構成例を説明する。まず、最適化装置10を含む最適化システムの例を説明する。
【0178】
図12は、最適化システムの例を示す図である。
最適化システム100は、最適化装置10および情報処理装置110を有する。
情報処理装置110は、CPU(Central Processing Unit)111、メモリ112およびバス113を有する。
【0179】
CPU111は、情報処理装置110を制御するプロセッサである。CPU111は、メモリ112に記憶されたプログラムを実行する。メモリ112は、例えば、RAM(Random Access Memory)であり、CPU111により実行されるプログラムやCPU111により処理されるデータを記憶する。バス113は、PCIe(Peripheral Component Interconnect Express)などの内部バスである。バス113には、CPU111、メモリ112および最適化装置10が接続される。
【0180】
CPU111は、バス113を介して、最適化装置10に対する温度Tの設定や規定回数C1,C2に関する情報の入力を行う。また、CPU111は、制約識別符号(識別コード、上限値、下限値を含む)や不等式制約毎の係数λを、最適化装置10に入力する。制御部18は、こうして入力された情報に基づいて、探索部D1に対する設定を行う。例えば、ユーザは、情報処理装置110に接続された入力部(図示を省略している)を操作して、制約識別符号や係数λを、情報処理装置110に設定することができる。
【0181】
このように、最適化装置10は、情報処理装置110に接続されて、不等式制約を含む組合せ最適化問題に対する演算をハードウェアにより高速に実行するアクセラレータとして利用される。
【0182】
次に、最適化装置の他の回路構成例を説明する。
図13は、最適化装置の他の回路構成例を示す図である。
最適化装置10aは、探索部D2および制御部18aを有する。探索部D2は、状態保持部11a、レジスタ11b1〜11bN、h計算部11c1〜11cN、ΔE計算部12a1〜12aN、ΔC加算部13d1〜13dN、判定部14a1〜14aNおよび選択部14bを有する。
【0183】
ここで、図13では、h計算部11c1〜11cN、ΔE計算部12a1〜12aNおよびΔC加算部13d1〜13dNについて、インデックスiの状態変数に対応することが分かり易いように「h」計算部などのように添え字iを付して名称を表記している。また、状態保持部11a、レジスタ11b1〜11bNおよびh計算部11c1〜11cNは、前述の状態保持部11に対応する。
【0184】
状態保持部11aは、状態ベクトルを保持する。状態保持部11aは、選択部14bから供給されるupdate信号に基づいて、状態ベクトルのうちの1つの状態変数の値を反転させることで、状態ベクトルを更新する。update信号は、反転させる状態変数のインデックスを示す。
【0185】
レジスタ11b1、h計算部11c1、ΔE計算部12a1、ΔC加算部13d1および判定部14a1が1番目の状態変数に関する演算を行う。また、レジスタ11b2、h計算部11c2、ΔE計算部12a2、ΔC加算部13d2および判定部14a2が2番目のスピンビットに関する演算を行う。同様に、「11b1」、「12a1」などの符号の末尾の数値iがインデックスiの状態変数に対応する演算を行うことを示す。すなわち、探索部D2は、レジスタ、h計算部、ΔE計算部、ΔC加算部および判定部のセット(1スピンビットに関する演算を行う演算処理回路の一単位であり、「ニューロン」と呼ばれることもある)を、N個有する。N個のセットが並列に、各セットに対応する状態変数に関する演算を行う。
【0186】
なお、あるニューロンにおけるΔC加算部には、他のニューロンにおけるレジスタに保持される重み値およびh計算部により保持される局所場が入力されるが、当該ΔC加算部と他のニューロンのレジスタおよびh計算部との信号線の図示を省略している。
【0187】
例えば、ΔC加算部13d1の例でいえば、W11,W21,…,WN1のうち、不等式制約用変数x(k∈Ω)に対応する1以上の重み値(Wk1=ak1)がΔC加算部13d1に供給される。また、h計算部11c1〜11cNのうち、不等式制約用変数xに対応する1以上の局所場(h)がΔC加算部13d1に供給される。なお、状態変数x(i∈Ψ)に対応するWi1やhがΔC加算部13d1に供給されても、図8で示した識別コードに基づく制御によりΔC加算部13d1は、ΔCi1=0とする。
【0188】
また、ΔC加算部13d2の例でいえば、W12,W22,…,WN2のうち、不等式制約用変数x(k∈Ω)に対応する1以上の重み値(Wk2=ak2)がΔC加算部13d2に供給される。また、h計算部11c1〜11cNのうち、不等式制約用変数xに対応する1以上の局所場(h)がΔC加算部13d2に供給される。なお、状態変数x(i∈Ψ)に対応するWi2やhがΔC加算部13d2に供給されても、識別コードに基づく制御によりΔC加算部13d2は、ΔCi2=0とする。
【0189】
以下では、主に、レジスタ11b1、h計算部11c1、ΔE計算部12a1、ΔC加算部13d1および判定部14a1を例示して説明する。同名の構成であるレジスタ11b2〜11bN、h計算部11c2〜11cN、ΔE計算部12a2〜12aN、ΔC加算部13d2〜13dNおよび判定部14a2〜14aNも同様の機能である。
【0190】
レジスタ11b1は、状態変数xとそれ以外の他の状態変数との間の重み値W1j(j=1〜N)を記憶する記憶部である。ここで、状態変数の数Nに対して、重み値の総数はN^2である。レジスタ11b1には、N個の重み値が格納される。
【0191】
レジスタ11b1は、インデックス1に対して、N個の重み値W11,W12,…,W1Nを記憶する。なお、Wii=W11=0である。レジスタ11a1は、選択部14bにより供給されるindex=jに対応する重み値W1jをh計算部11c1に出力する。
【0192】
h計算部11c1は、レジスタ11b1から供給される重み係数W1jを用いて、式(3)、(4)に基づく局所場hを計算する。h計算部11c1は、前回計算された局所場hを保持するレジスタを有し、index=jで示される状態変数の反転方向に応じたδh(j)を、hに積算することで、当該レジスタに格納されるhを更新する。hの初期値は、問題に応じて、h計算部11c1のレジスタに予め設定される。また、bの値も、問題に応じて、h計算部11c1のレジスタに予め設定される。h計算部11c1は、計算した局所場hをΔE計算部12a1に出力する。
【0193】
ΔE計算部12a1は、h計算部11c1から供給される局所場hに基づいて、式(2)によりΔEを計算し、ΔC加算部13d1に出力する。ここで、ΔE計算部12a1は、状態変数xの反転方向に応じてhの符号を決めることでΔEを計算する。
【0194】
ΔC加算部13d1は、式(15)に基づいて不等式制約に対する超過量の変化値に応じたペナルティ値ΔCを計算し、ΔE計算部12a1から供給されるΔEに対して、ΔCを加算する。ΔC加算部13d1は、ΔE+ΔCを判定部14a1に出力する。
【0195】
判定部14a1は、ΔC加算部13d1から供給されるΔE+ΔCに基づいて、反転可否を示すフラグFを選択部14bに出力する。
選択部14bは、判定部14a1〜14aNから供給されるフラグF〜Fに基づいて、反転対象の状態変数のインデックスjを選択し、状態保持部11aおよびレジスタ11b1〜11bNに供給する。なお、選択部14bは、判定部14a1〜14aNに対して、反転対象の状態変数のインデックスjを選択したか否かを示すフラグを供給することで、オフセット値の生成を促すが、当該フラグを供給する信号線の図示を省略している。
【0196】
制御部18aは、探索部D2における判定部14a1〜14aNに対する温度Tの設定、状態変数の更新回数の制御、ΔE計算部12a1〜12aNおよびΔC加算部13d1〜13dNに対する制約識別符号の入力を行う。ΔE計算部12a1〜12aNのΔE計算およびΔC加算部13d1〜13dNのΔC計算は、図8で示されるように、制約識別符号により制御される。
【0197】
最適化装置10aも、最適化装置10と同様に、不等式制約を含む最適化問題の求解を効率的に実行できる。最適化装置10aを用いることで、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。すなわち、最適化装置10aの応用範囲を広げることができ、より多くの問題の求解に、最適化装置10aを利用可能になる。
【0198】
[第2の実施の形態]
次に、第2の実施の形態を説明する。前述の第1の実施の形態と相違する事項を主に説明し、共通する事項の説明を省略する。
【0199】
第1の実施の形態では、最適化装置10または最適化装置10aが単一の探索部を有する構成を例示した。一方、複数の探索部を有する最適化装置において、レプリカ交換法を用いて、最適化問題の求解を行う構成も考えられる。
【0200】
図14は、第2の実施の形態の最適化装置の例を示す図である。
最適化装置50は、探索部51a1〜51aM(M(Mは2以上の整数)個の探索部)および交換制御部52を有する。
【0201】
探索部51a1〜51aMは、第1の実施の形態の探索部D1または探索部D2と同様の回路構成を有する。図14では、探索部51ajの例が示されている。
探索部51ajは、状態保持部11、エネルギー変化計算部12、ペナルティ加算部13、更新制御部14、制約超過量計算部15および制約識別符号入力部16を有する。状態保持部11、エネルギー変化計算部12、ペナルティ加算部13、更新制御部14、制約超過量計算部15および制約識別符号入力部16のそれぞれの機能は、第1の実施の形態における同名の要素と同様である。また、探索部51ajは、エネルギー計算部17を更に有する(図示を省略している)。
【0202】
交換制御部52は、探索部51a1〜51aMのそれぞれに互いに異なる温度値を設定し、探索部51a1〜51aMを用いたレプリカ交換法による基底状態探索を制御する。交換制御部52は、探索部51a1〜51aMのそれぞれによる基底状態探索の繰り返し回数到達または一定時間経過後に、探索部間で温度値または各探索部に保持される状態変数の値を入れ替える。例えば、交換制御部52は、温度が隣接する2つの探索部のペアで、所定の確率(交換確率)により温度を交換することが考えられる。交換確率としては、例えば、メトロポリス法に基づく確率が用いられる。あるいは、交換制御部52は、温度の代わりに、温度が隣接する2つの探索部のペアで、当該2つの探索部が保持する状態(複数の状態変数の値)や各状態変数に対応する局所場を交換してもよい。
【0203】
このように、最適化装置50によれば、レプリカ交換法を用いる場合でも、不等式制約を含む最適化問題の求解を効率的に実行できる。最適化装置50によれば、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。最適化装置50の応用範囲を広げることができ、より多くの問題の求解に、最適化装置50を利用可能になる。
【符号の説明】
【0204】
10 最適化装置
11 状態保持部
12 エネルギー変化計算部
13 ペナルティ加算部
14 更新制御部
14a 遷移可否判定部
14b 選択部
15 制約超過量計算部
16 制約識別符号入力部
図1
図2
図3
図4
図5
図6
図7
図8
図9
図10
図11
図12
図13
図14