JP2021184148A - 最適化装置、最適化方法、および最適化プログラム - Google Patents

最適化装置、最適化方法、および最適化プログラム Download PDF

Info

Publication number
JP2021184148A
JP2021184148A JP2020088882A JP2020088882A JP2021184148A JP 2021184148 A JP2021184148 A JP 2021184148A JP 2020088882 A JP2020088882 A JP 2020088882A JP 2020088882 A JP2020088882 A JP 2020088882A JP 2021184148 A JP2021184148 A JP 2021184148A
Authority
JP
Japan
Prior art keywords
value
replica
state
replicas
interaction
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
JP2020088882A
Other languages
English (en)
Inventor
悟 半田
Satoru Handa
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujitsu Ltd
Original Assignee
Fujitsu Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujitsu Ltd filed Critical Fujitsu Ltd
Priority to JP2020088882A priority Critical patent/JP2021184148A/ja
Priority to EP21153978.8A priority patent/EP3913506A1/en
Priority to US17/168,557 priority patent/US11790130B2/en
Priority to CN202110183741.8A priority patent/CN113705851A/zh
Publication of JP2021184148A publication Critical patent/JP2021184148A/ja
Withdrawn legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N5/00—Computing arrangements using knowledge-based models
    • G06N5/01—Dynamic search techniques; Heuristics; Dynamic trees; Branch-and-bound
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/08—Computing arrangements based on specific mathematical models using chaos models or non-linear system models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/044—Recurrent networks, e.g. Hopfield networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/047—Probabilistic or stochastic networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Databases & Information Systems (AREA)
  • Business, Economics & Management (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Computer Hardware Design (AREA)
  • Nonlinear Science (AREA)
  • Geometry (AREA)
  • Development Economics (AREA)
  • Marketing (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Game Theory and Decision Science (AREA)
  • Computational Linguistics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

【課題】複数のレプリカを用いた場合の解探索能力を向上させる。【解決手段】最適化装置は、複数の状態変数を有する複数のレプリカそれぞれについて、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内でのそのレプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定する。そして最適化装置は、第１状態変数の値を更新した場合における相互作用の強さの変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、第１状態変数の値を更新するか否かを決定する。【選択図】図１８

Description

本発明は、最適化装置、最適化方法、および最適化プログラムに関する。
ノイマン型コンピュータが不得意とする問題として、大規模な離散最適化問題がある。離散最適化問題を計算する装置としては、例えば、イジング型の評価関数（エネルギー関数などとも呼ばれる）を用いたイジングマシン（ボルツマンマシンとも呼ばれる）がある。
イジングマシンによる計算では、計算対象の問題は磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えられる。そして、マルコフ連鎖モンテカルロ法により、イジング型の評価関数の値（イジングモデルのエネルギーに相当する）が最小となる状態の探索が行われる。以下、マルコフ連鎖モンテカルロ法を、ＭＣＭＣ（Markov-Chain Monte Carlo）法と略す。ＭＣＭＣ法では、例えばメトロポリス法またはギブス法で規定される状態遷移の受け入れ確率で、その状態遷移が受け入れられる。
ＭＣＭＣ法の一種として、レプリカ交換法（交換モンテカルロ法またはパラレルテンパリング（parallel tempering）法とも呼ばれる）がある。レプリカ交換法は複数の温度を用いたＭＣＭＣ処理を互いに独立に行い、ある試行回数ごとに、各ＭＣＭＣ処理で得られるエネルギーを比較し、適切な確率で２つの温度に対する状態を交換するという操作を行う方法である。レプリカ交換によれば、温度を徐々に下げていく疑似焼き鈍し法と比べて、局所解に拘束される可能性が抑えられ、全探索空間を効率よく探索できる。
解探索に関する技術としては、さらに、探索過程における状態のエネルギー分布をボルツマン分布に近づけることにより、遷移候補探索の並列化による収束特性の変化、解の精度の悪化を抑制する最適化装置も提案されている。レプリカ交換に関する技術としては、例えば、計算時間を短縮するイジングマシンが提案されている。またレプリカ交換法による擬似焼鈍し動作を行う最適化装置の回路規模の増大を抑制する最適化装置が提案されている。さらに回路の物量を削減しつつ、メトロポリス法に基づく確率的な処理を可能とする情報処理装置も提案されている。複数のレプリカを用いて解探索を行う技術としては、Collective Monte Carlo（ＣＭＣ）と呼ばれる方法や、Robust Ensemble（ＲＥ）と呼ばれる方法も提案されている。
特開２０２０−０６１０１２号公報 特開２０１８−００５５４１号公報 特開２０１９−１９４７６５号公報 特開２０１９−０８２７９３号公報
Gregoire Clarte and Antoine Diez,"Collective sampling through a Metropolis-Hastings like method: kinetic theory and numerical experiments", arXiv:1909.08988v1 [math.ST], 18 Sep. 2019 Baldassi, Carlo. et. al., "Unreasonable Effectiveness of Learning Neural Networks: From Accessible States and Robust Ensembles to Basic Algorithmic Schemes", PNAS E7655-E7662, Published online 15 Nov. 2016
ＭＣＭＣ法を高速化するために、多数のレプリカを用いて集団で探索をする様々な方法が提案されているが、いずれの方法においても、集団探索の効果が十分には発揮できない場合がある。例えば、遷移先の候補の選び方が１−ｂｉｔフリップ（複数のビットのうちの１つの値を反転させる）の場合、各ビットが反転対象として等確率で選択され、選択されたビットを反転させた状態への遷移確率は遷移前後のエネルギー差のみで決定される。そのため、どのレプリカもエネルギー勾配のみに従って状態変化し、状態遷移の過程が同じ道筋を辿ってしまう可能性がある。その結果、複数のレプリカが同じ局所解に留まってしまい、十分に広く状態空間を探索できない場合がある。
なお、このような問題は、状態変数の値が離散的な場合に限らず、状態変数が連続的な値を取り得る最適化問題においても同様に生じる。
１つの側面では、本発明は、複数のレプリカを用いた場合の解探索能力を向上させることを目的とする。
１つの案では、以下に示す記憶部と処理部とを有する最適化装置が提供される。
記憶部は、複数のレプリカそれぞれの複数の状態変数の値を記憶する。処理部は、複数のレプリカそれぞれについて、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内でのレプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定する。そして処理部は、第１状態変数の値を更新した場合における相互作用の強さの変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、第１状態変数の値を更新するか否かを確率的に決定する。
１態様によれば、複数のレプリカを用いた場合の解探索能力を向上させることができる。
第１の実施の形態に係る最適化方法の一例を示す図である。 第２の実施の形態のシステム構成の一例を示す図である。 サーバのハードウェアの一例を示す図である。 イジングマシンの一例を示す図である。 イジングモデルの模式図である。 レプリカ交換の一例を示す図である。 Ｏｎｅ−Ｈｏｔ制約下での１−ｂｉｔフリップの例を示す図である。 ２Ｗ１Ｈの制約を説明する図である。 イジングマシンの解探索機能の一例を示す図である。 解探索エンジンにおける処理の一例を示す図である。 解探索処理の手順の一例を示すフローチャートである。 レプリカごとの解探索処理の手順の一例を示すフローチャートである。 レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。 第１の更新ビット選択方法による更新ビット選択処理の手順の一例を示すフローチャートである。 第２の更新ビット選択方法の処理手順の一例を示すフローチャートである。 更新ビット選択のためのツリー状に接続されたセレクタの一例を示す図である。 第３の更新ビット選択方法の処理手順の一例を示すフローチャートである。 レプリカ間に斥力の相互作用を設定した場合のエネルギーランドスケープを示す図である。 レプリカ間に引力の相互作用を設定した場合のエネルギーランドスケープを示す図である。 第１の検証例を示す図である。 第２の検証例を示す図である。
以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
図１は、第１の実施の形態に係る最適化方法の一例を示す図である。図１には、解探索方法を実現する最適化装置１０が示されている。最適化装置１０は、ノイマン型コンピュータであってもよく、非ノイマン型コンピュータであってもよい。例えば最適化装置１０は、最適化用のプログラムを実行することにより、最適化方法を実施することができる。また最適化装置１０は、イジングモデルを用いた最適化問題の求解を行うイジングマシンであってもよい。イジングマシンには、量子ビットを用いた量子コンピュータ、量子ビットの量子現象をデジタル回路上で再現した装置などが含まれる。
最適化装置１０は、記憶部１１と処理部１２とを有する。記憶部１１は、例えば最適化装置１０が有するメモリ、またはストレージ装置である。処理部１２は、例えば最適化装置１０が有するプロセッサ、または演算回路である。演算回路には、量子ビットまたは量子ビットの仕組みを再現するニューロン回路が含まれる。
記憶部１１は、複数のレプリカ２〜４それぞれの複数の状態変数の値を記憶する。
処理部１２は、複数のレプリカ２〜４を用いて、最適化問題を求解する。例えば処理部１２は、最適化問題に応じて定義された目的関数の値が最小となる状態変数の値を求める。目的関数は、最適化問題を表すモデルのエネルギーと呼ばれることもある。最適化問題がイジングモデルで表される場合、イジングモデルのハミルトニアンが、エネルギーを示す目的関数に相当する。
解探索のために処理部１２は、複数のレプリカ２〜４それぞれについて状態遷移（状態変数の値を更新すること）を繰り返し、生成された状態における複数の状態変数の値に基づいて目的関数の値を計算する。その際、処理部１２は、レプリカ間の相互作用を考慮して、レプリカの状態遷移を行う。レプリカ間の相互作用としては、例えばレプリカ間の距離に応じた引力または斥力が考えられる。ｋ番目のレプリカｘkとｌ番目のレプリカｘlとの距離は、ｄ（ｘk，ｘl）と表記する（ｋ，ｌは、１以上の整数）。例えば処理部１２は、複数のレプリカ２〜４それぞれについて以下のようにして状態遷移を行う。
処理部１２は、レプリカが有する複数の状態変数のうちの第１状態変数の値を更新した場合における、複数の状態変数の値の組み合わせが取り得る空間を示す状態空間１内でのそのレプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定する。相互作用の強さは、例えば他のレプリカそれぞれとの距離の合計に基づく値である。相互作用の強さは、相互作用のエネルギーＧ（Ｘ）と呼ぶこともできる。相互作用の強さは、例えば後述の式（１５）または式（１６）で表すことができる。ｌ番目のレプリカのｊ0番目の状態変数を更新した場合の相互作用の強さの変化量は、ΔＧ＝Ｇ（ｘl［ｊ0］）−Ｇ（ｘl）と表すことができる。
そして処理部１２は、第１状態変数（例えばｊ0番目の状態変数）の値を更新した場合における、第１状態変数の値を更新するか否かを決定する。この決定は、相互作用の強さの変化量に応じた提案確率（ｇ（ｘl→ｘl［ｊ0］）と、目的の確率分布に応じた受け入れ確率（ａ（ｘl→ｘl［ｊ0］）とに基づいて、確率的に行われる。目的の確率分布は、例えばギブス分布である。提案確率と受け入れ確率とに基づいてレプリカが状態遷移する遷移確率は、例えばメトロポリスヘイスティング法に従う。
処理部１２は、第１状態変数の値を更新すると決定した場合、第１状態変数の値を更新した後のレプリカの複数の状態変数の値に基づいて、目的関数の値を計算する。また処理部１２は、記憶部１１内におけるレプリカの第１状態変数の値を更新する。そして処理部１２は、複数のレプリカ２〜４それぞれの複数の状態変数の一の状態変数の値の更新を繰り返し、目的関数の値が所定の条件を満たしたときの複数の状態変数の値を出力する。例えば処理部１２は、複数のレプリカ２〜４の更新を所定回数繰り返した後、目的関数の値が最小となる複数の状態変数の値の組み合わせを出力する。
このようにして、レプリカ間の相互作用を考慮したレプリカの状態遷移により、解の探索が行われる。すなわち最適化装置１０は、レプリカ間の相互作用を考慮したことで、複数のレプリカ２〜４による状態空間１内を網羅的に探索することができる。しかもメトロポリスヘイスティング法を用いることで、最適化装置１０は、レプリカ間の相互作用による影響を適切な形で計算に組み込むことができる。
なお処理部１２は、状態空間１に対して適切な距離を定義し、レプリカ間の距離を定める。そして処理部１２は、その距離を用いてレプリカ間の相互作用の強さを決定し、メトロポリスヘイスティング法における遷移候補先の分布（提案分布）を定義し、計算に組み込む。メトロポリスヘイスティング法は提案分布が非対称な場合に対応している。そのため、提案分布の決め方に自由度がある。そこで、処理部１２は、メトロポリスヘイスティング法における提案分布（提案確率の定義）の自由度を利用し、提案確率内にレプリカ間相互作用を導入している。
レプリカ間の相互作用として、例えば斥力の相互作用を生じさせることができる。この場合、処理部１２は、第１状態変数の値を更新すると状態遷移判断対象のレプリカと他のレプリカとの距離が遠ざかる場合に、相互作用の強さを増加させる。処理部１２は、相互作用の強さの増加量が大きいほど提案確率を大きくする。そして処理部１２は、提案確率が大きい状態変数ほど、値を更新する候補として選択される確率を高くする。その結果、例えば複数のレプリカ２〜４が、同じ局所解に嵌まり出られなくなることを抑止できる。
またレプリカ間の相互作用として、引力の相互作用を生じさせることもできる。この場合、処理部１２は、第１状態変数の値を更新すると状態遷移判断対象のレプリカと他のレプリカとの距離が近づく場合に、相互作用の強さを増加させる。処理部１２は、相互作用の強さの増加量が大きいほど提案確率を大きくする。これにより、局所解に嵌まって出られなくなっているレプリカを、他のレプリカからの引力により、局所解から抜け出させることが可能となる。
状態空間１が離散的であり、状態変数の値が二値（例えば「１」または「０」）のみを取り得る場合、２つのレプリカ間の距離として、例えばハミング距離（またはその単調増加関数）を用いることができる。この場合、処理部１２は、すべてのレプリカ間のハミング距離を定義し、それによりレプリカ相互作用の強さを計算する。レプリカ間の距離は、後述の式（１９）のように表すことができる。
なお処理部１２は、第１状態変数の値を更新することについての提案確率を、例えば規格化定数によって規格化する。例えば第１状態変数の値を更新した場合における相互作用の強さの変化量をΔＧ、レプリカに設定された温度パラメータの値の逆数である逆温度をβとする。このとき処理部１２は、ｅｘｐ（−βΔＧ）を所定の規格化定数で除算した値を提案確率とする。この提案確率は、例えば後述の式（１７）で表すことができる。ｅｘｐ（−βΔＧ）はギブス分布を表しており、提案確率の定義にギブス分布を用いたことで、目的関数（エネルギー）においてギブス分布を保つことが容易になる。
また処理部１２は、１とｅｘｐ（−βΔＧ）とのうちの小さい方を、所定の規格化定数で除算した値を提案確率としてもよい。この提案確率は、例えば後述の式（１８）で表すことができる。これによりｅｘｐ（−βΔＧ）が１を超える場合は１とみなされることとなり、相互作用の強さの変化量が状態変数間で大きく異なる場合の提案確率への影響の差を弱めることができる。
ここで、規格化定数について説明する。従来の提案分布は複数の状態変数が等確率（１／Ｎ）で遷移候補として選択される（Ｎは、状態変数の数を示す１以上の整数）。この場合、規格化定数はＮ（それぞれの遷移先の重みが共通で１ということ）である。図１の最適化装置１０では、遷移候補となる状態変数それぞれの遷移確率が異なり、規格化定数が遷移前の現在状態に依存する。そこで処理部１２において規格化定数を計算することとなる。
例えば処理部１２は、複数の状態変数それぞれを第１状態変数とした場合における複数の状態変数それぞれについてのｅｘｐ（−βΔＧ）の値の総和を、規格化定数とする。この規格化定数は、例えば後述の式（２３）で表すことができる。なお、相互作用がハミング距離の一次関数の場合、処理部１２は、レプリカの状態遷移ごとに状態遷移前と状態遷移後との規格化定数の差分計算を行い、差分の累積値を計算し（累積計算）、最新の規格化定数とすることができる。ハミング距離の一次関数は、後述の式（１９）に示すような関数である。
規格化定数の累積計算をする場合、処理部１２は、レプリカを状態遷移させるごとに、更新させる状態変数の決定に使用した規格化定数を記憶部１１に格納しておく。そして処理部１２は、レプリカの状態遷移の際に使用した規格化定数の値と、前回の状態遷移の前後で生じる規格化定数の値の差分とに基づいて、今回の状態遷移で使用する規格化定数の値を算出する。前回の状態遷移の前後で生じる規格化定数の値の差分は、例えば後述の式（２４）で表される。これにより、規格化定数を効率的に算出することができる。
なお処理部１２は、他のレプリカそれぞれとの距離の平方根の合計に基づく値を、相互作用の強さとしてもよい。この場合の相互作用の強さは、例えば後述の式（１６）で表される。これにより、距離が遠い他のレプリカよりも距離が近い他のレプリカからの相互作用が、相対的により強く働くようにすることができる。例えば、複数のレプリカ２〜４が同じ局所解に嵌まることを抑止する場合、その局所解の近辺に存在するレプリカ間に強い斥力を働かせることで、局所解からの脱出を促進することができる。この場合、その局所解とは遠く離れた位置のレプリカからの影響が少ない方が、局所解からの脱出が容易となる。
また処理部１２は、複数の状態変数の中から、値の更新を受け入れることができる状態変数を先に特定し、その中から今回のレプリカの状態遷移において値を更新する状態変数を決定することもできる。この場合、処理部１２は、複数の状態変数それぞれについて、受け入れ確率に基づいて、状態変数の更新が提案された場合に更新を受け入れるか否かを確率的に決定する。そして処理部１２は、更新を受け入れると判定された状態変数の中から、提案確率が高い状態変数ほど選択される可能性を高くして、少なくとも１つの状態変数を更新対象に決定する。これにより、選択した状態変数の値の更新の棄却（更新を受け入れないとの判定）が繰り返され、値を更新する状態変数の決定に時間がかかることを抑止できる。
〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、目的関数の値が最小となる各状態変数の値の組み合わせを計算するイジングマシンを用いたシステムの例である。なお第２の実施の形態におけるイジングマシンは、第１の実施の形態に示した最適化装置１０の一例である。イジングマシンでは、求解対象の問題をイジングモデルで表し、そのイジングモデルのエネルギーが最小値となるビットの値の組み合わせを探索する。イジングモデルのエネルギーを計算する式（ハミルトニアン）が、目的関数である。
図２は、第２の実施の形態のシステム構成の一例を示す図である。サーバ１００には、ネットワーク２０を介して端末装置３１，３２，・・・が接続されている。端末装置３１，３２，・・・は、組み合わせ最適化問題の求解を依頼するユーザが使用するコンピュータである。サーバ１００は、端末装置３１，３２，・・・から組み合わせ最適化問題の求解の依頼を受け付け、組み合わせ最適化問題に対応するイジングモデルのエネルギー関数であるハミルトニアンを生成する。サーバ１００には、イジングマシン３００の制御装置２００が接続されている。サーバ１００は、生成したハミルトニアンを用いてエネルギーの最小値の探索要求を制御装置２００に入力する。
制御装置２００は、イジングマシン３００を制御し、サーバ１００から入力された探索要求に応じて、エネルギーの最小値の解探索を行う。例えば制御装置２００は、各ニューロンについての結合先のニューロンのｉｄを、結合先情報としてイジングマシン３００に送信する。また、制御装置２００は、ローカルフィールドの初期値（例えばバイアス係数）や、値が０ではない重み係数、アニーリング条件などについてもイジングマシン３００に送信する。
イジングマシン３００は、制御装置２００からの制御に基づいて、デジタル回路を用いたイジングモデルの状態遷移のシミュレーションを行い、エネルギーの最小値を探索する。
図３は、サーバのハードウェアの一例を示す図である。サーバ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。
メモリ１０２は、サーバ１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に利用する各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。
バス１０９に接続されている周辺機器としては、ストレージ装置１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７およびネットワークインタフェース１０８がある。
ストレージ装置１０３は、内蔵した記録媒体に対して、電気的または磁気的にデータの書き込みおよび読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、および各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。
グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、有機ＥＬ（Electro Luminescence）を用いた表示装置や液晶表示装置などがある。
入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。
光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取り、または光ディスク２４へのデータの書き込みを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。
機器接続インタフェース１０７は、サーバ１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。
ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。ネットワークインタフェース１０８は、例えばスイッチやルータなどの有線通信装置にケーブルで接続される有線通信インタフェースである。またネットワークインタフェース１０８は、基地局やアクセスポイントなどの無線通信装置に電波によって通信接続される無線通信インタフェースであってもよい。
サーバ１００は、以上のようなハードウェアによって、第２の実施の形態の処理機能を実現することができる。なお、制御装置２００も、サーバ１００と同様のハードウェアにより実現することができる。
図４は、イジングマシンの一例を示す図である。イジングマシン３００は、ニューロン回路３１１，３１２，…，３１ｎ、制御回路３２０、およびメモリ３３０を有する。
ニューロン回路３１１〜３１ｎはそれぞれ、自身以外の複数の他のニューロン回路との接続の有無を示す複数の重み係数の値と、複数の他のニューロン回路の複数の出力信号との積の総和に基づく第１の値を算出する。そしてニューロン回路３１１〜３１ｎそれぞれは、第１の値にノイズ値を加算した第２の値と閾値との比較結果に基づき、０または１のビット値を出力する。複数のレプリカを用いた解探索を行う場合、複数のニューロン回路を用いて１つのレプリカの解探索が行われる。
制御回路３２０は、制御装置２００から供給される情報に基づいて、イジングマシン３００の初期設定処理などを行う。また、制御回路３２０は、レプリカ交換を行う場合、２つのレプリカ間の温度パラメータの値の交換の有無を判断し、温度パラメータの値を交換する場合、各レプリカの解探索を行うニューロン回路に入力する温度パラメータの値を更新する。
さらに、制御回路３２０は、更新対象ニューロンを決定する処理が所定回数だけ繰り返された後、メモリ３３０に保持されている１レプリカの状態変数に対応する各ニューロンのビット値を取得し、最適化問題に対する解として制御装置２００に送信する。
制御回路３２０は、例えばＡＳＩＣやＦＰＧＡなどの特定用途の電子回路にて実現できる。なお、制御回路３２０は、ＣＰＵやＤＳＰなどのプロセッサであってもよい。その場合、プロセッサは、図示しないメモリに記憶されたプログラムを実行することで、上記の処理を行う。
メモリ３３０は、例えば各ニューロンのビット値を保持する。メモリ３３０は、例えばレジスタやＲＡＭなどによって実現できる。メモリ３３０には、エネルギーの最小値や最小値が得られたときの各ニューロンのビット値を保持しておくこともできる。この場合、制御回路３２０は、更新対象ニューロンを決定する処理が所定回数、繰り返されたのちに、エネルギーの最小値や最小値が得られたときの各ニューロンのビット値をメモリ３３０から取得して、制御装置２００に送信してもよい。
なお第１の実施の形態に示した最適化装置１０も、図４に示したイジングマシン３００と同様のハードウェアにより実現することができる。
次に、求解の対象となるイジング型の最小値求解問題（イジング型問題）について説明する。イジング型問題はイジングモデルで表される。
図５は、イジングモデルの模式図である。イジングモデル３０は、格子状に複数のビット３１が配置される。ビット３１は、磁石を模したものであり、スピンとも呼ばれる。隣接するビット間には相互作用が働く。相互作用の大きさは、重み係数で表される。イジングモデル３０のエネルギーは、次の式（１）で表される。
Figure 2021184148
右辺の１項目は、Ｎ個の状態変数の全組み合わせについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘiはｉ番目の状態変数、ｘjはｊ番目の状態変数を表し、Ｗijは、ｘiとｘjとの結合の強さを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂi）の総和を求めたものである。Ｗijが正の場合、ｘiとｘjとが同じ値となるような相互作用が働く。またＷijが負の場合、ｘiとｘjとが異なる値となるような相互作用が働く。なおＷij＝Ｗjiであり、Ｗii＝０である。
最小値求解問題は、式（１）で与えられるエネルギーの最小値を求める問題である。イジングマシン３００は、このような最小値給回問題についてＭＣＭＣを用いて解く。例えばイジングマシン３００は、ビットを１つ反転した場合のエネルギー変化を計算する。ｉ番目のビットを反転させた場合「ｘi→ｘi´（δｘi＝ｘi´−ｘi）」、エネルギーの変化値は式（２）で表される。
Figure 2021184148
式（２）の右辺の括弧内の式は、ｉ番目のビットの局所場（総入力）を表している。出力変化分δｘiと局所場の符号が一致すればエネルギーは減少する。イジングマシン３００は、エネルギーの変化値ΔＥiの増減に応じて、ｉ番目のビットの反転を受け入れるか否かを決定する。なお、式（２）は、１ビットのみを反転させる場合にのみ正しい式である。
エネルギー増分を示す式（２）は、以下のように書き換えることができる。
Figure 2021184148
Figure 2021184148
ｈiはｉ番目のビットの局所場である。ｊ番目のビットｘjが反転したときのｉ番目のビットの局所場ｈiの変化分δｈi (j)は、以下の式（５）で表される。
Figure 2021184148
局所場ｈiを保存するレジスタを用意し、ｊ番目のビットが反転したときに、式（５）に示す値を、保存しておいた局所場ｈiに加算することで、常に正しいｈiが得られる。
以上のような計算により、ｉ番目のビットの反転させた場合のエネルギーの増分を求めることができる。イジングマシン３００は、得られたエネルギーの増分に基づいて、ｉ番目のビットの反転を受け入れるか否かを判断する。例えばイジングマシン３００は、メトロポリス法に従って、ビットの反転を受け入れるか否かを判断する。メトロポリス法に従う場合、エネルギーの増分が負（エネルギーが減少する）であれば、ビットの反転が受け入れられる。エネルギーの増分が正（エネルギーが増加する）であれば、エネルギーの増分に応じた確率によって、ビットの反転が受け入れられるか否かが判断される。
エネルギーの増分が正の場合にビットの反転が受け入れられる確率は、温度パラメータを用いて調整することができる。例えばイジングマシン３００は、温度パラメータの値が大きいほど、エネルギーの増分が正の場合にビットの反転を受け入れる確率を高くする。これにより、温度パラメータの値を大きくすることで、イジングモデルのエネルギーの状態が局所解から抜け出す可能性を高くすることができる。
温度パラメータをＴとしたとき、逆温度β＝１／Ｔをとする。例えばイジングマシン３００は、エネルギーの変化値ΔＥijと逆温度βを用いて、ｉ番目の状態変数の状態遷移の受け入れ確率を以下の式（６）により決定することで、確率的探索を行う。
Figure 2021184148
式（６）における関数ｆ（ｘ）は、メトロポリス法では以下の式（７）となる。
Figure 2021184148
なお温度パラメータの値が大きいと、局所的な探索が難しくなる。そこでイジングマシン３００は、例えば温度パラメータの値が異なる複数のレプリカを用いて解探索を行う。この場合、イジングマシン３００はレプリカ交換を行ってもよい。
図６は、レプリカ交換の一例を示す図である。レプリカ交換では、複数のレプリカが用いられる。レプリカは、求解対象の問題の状態変数の集合のコピーである。イジングマシン３００は、それぞれのレプリカの温度パラメータに異なる値を設定する。図６の例では、４つのレプリカそれぞれに、Ｔ1，Ｔ2，Ｔ3，Ｔ4の温度パラメータが設定されている（Ｔ1＜Ｔ2＜Ｔ3＜Ｔ4）。
イジングマシン３００は、複数のレプリカそれぞれについて、ＭＣＭＣにより状態変化させる。そしてイジングマシン３００は、温度パラメータの値で並べたときに隣り合うレプリカ間で、所定の確率に従って温度パラメータの値を交換する。すると、各レプリカは温度軸方向にランダムウォークする。レプリカがランダムウォークをすることで、局所解に嵌まっても高温側に移動したときに局所解から脱出できる可能性がある。またレプリカが低温側に移動すると、局所的な探索を行うこともできる。
レプリカ交換のように多くのレプリカを用いて集団探索を行うことで、モンテカルロ法による解探索を高速化することができる。しかし単に複数のレプリカを用いて集団探索をしただけでは、複数のレプリカが同じ局所解に留まってしまい、十分に広く状態空間を探索できないという問題を解決できない。例えばイジングモデルの状態変数（ビット）の数がＮのとき、状態空間には２N個の状態が存在する。そのため状態変数の数が多くなると、実用上可能な数のレプリカで探査しても、集団探索の利益を享受するのは困難である。
そこでイジングマシン３００では、レプリカ間の距離に応じた相互作用を利用して、レプリカの状態を遷移させることで、状態空間内の効率的な探索を行う。これにより、複数のレプリカを用いた集団探索による解探索性能が向上する。
例えばレプリカ交換を行えば状態空間内を広範囲に探索できるが、レプリカ間の相互作用を考慮しない場合、各レプリカは、そのときの温度パラメータの値に応じて、独立にビットフリップ（マルコフ連鎖）を行うだけである。レプリカ間の相互作用を利用すれば、個々のレプリカのマルコフ連鎖において、複数のレプリカが同じ局所解に同時に留まることを抑止可能である。
また１−ｂｉｔフリップの場合、遷移先の候補の選び方としてＮ個のビットを等確率で選択すると、遷移確率はエネルギーの変化値ΔＥiのみで決定される。この場合、どのレプリカもエネルギー勾配のみに従って状態変化するため、同じ道筋を辿ってしまい、十分に広く状態空間を探索できていない可能性が大きい。さらに、どのレプリカも同じ局所解（すべてのビットｉでΔＥi＞０）にはまってしまったときに脱出するのも困難である。
イジングマシン３００では、遷移を受け入れるか否かの受容確率の計算にメトロポリス法ではなく、メトロポリス・ヘイスティングス法を用いる。これにより、レプリカ間の相互作用による影響を適切な形で計算に組み込むことができる。
例えば現在の状態ｘから次の状態Ｘ′を提案する確率をｇ（Ｘ→Ｘ′）とし、この状態遷移が受け入れられる確率をＡ（Ｘ→Ｘ′）とする。状態Ｘから状態Ｘ′に遷移する確率Ｗ（Ｘ→Ｘ′）は、以下の式（８）で得られる。
Figure 2021184148
目的の確率分布（例えばギブス分布）を表す関数（目的関数）をπ（Ｘ）とすると詳細つり合いの条件は以下の通りとなる。
Figure 2021184148
Figure 2021184148
式（１０）から、詳細つり合いを満たす受け入れ確率は、式（１１）の通りとなる。
Figure 2021184148
メトロポリス・ヘイスティングス法を適用した場合、受け入れ確率は以下の式（１２）で与えられる。
Figure 2021184148
この受け入れ確率は、提案確率が非対称でｇ（Ｘ→Ｘ′）≠ｇ（Ｘ′→Ｘ）の場合でも詳細つり合いの条件は満たされる。また提案確率が対称でｇ（Ｘ→Ｘ′）＝ｇ（Ｘ′→Ｘ）の場合は式（１３）のようなメトロポリスの受け入れ確率が得られる。
Figure 2021184148
ここで１−ｂｉｔフリップを考えるとき、レプリカ間の相互作用を考慮しない場合、Ｎ個のビットが等確率で反転の候補として選択され、提案確率は式（１４）となる。
Figure 2021184148
なお、メトロポリス・ヘイスティングス法は提案確率で示される提案分布が非対称な場合に対応している。そのため、提案分布の決め方に自由度がある。そこでイジングマシン３００は、提案確率にレプリカ間相互作用を導入する。
例えばイジングマシン３００は、離散空間である状態空間に対して適切な距離を定義し、レプリカ間の距離を定める。イジングマシン３００は、レプリカ間距離を用いてレプリカ相互作用を決め、メトロポリス・ヘイスティングス法における遷移候補先の分布（提案分布）を定義し、受け入れ確率の計算に組み込む。
レプリカ間距離の一例として、２つのレプリカの状態のハミング距離（或いはその単調増加関数）がある。イジングマシン３００は、すべてのレプリカ間のハミング距離を定義し、それによりレプリカ相互作用を導入する。
式（１４）に示すような提案分布は１／Ｎの等確率で遷移候補が選択されていたため、規格化定数はＮ（それぞれの遷移先の重みが共通で１ということ）である。レプリカ間の相互作用を導入した場合は遷移候補の重みが異なり、規格化定数が遷移前の現在状態に依存する。イジングマシン３００は、規格化定数も計算しなければならないが、相互作用がハミング距離の一次式の場合は、差分計算（累積計算）によって容易に計算できる。
以下、レプリカ間距離を考慮した提案確率の計算方法について具体的に説明する。まず提案確率の一般系を以下のように定義する。
Ｍ個（Ｍは１以上の整数）のレプリカからなる系を考える。１番目のレプリカの状態変数をｘl＝（ｘ1 l，ｘ2 l，・・・，ｘN l），ｘj l∈｛０，１｝とする。２つのレプリカｘlとｘkの距離（の増加関数）をｄ（ｘl，ｘk）とし、相互作用のエネルギーをＧ（ｘ）と与える。相互作用のエネルギーは、例えば式（１５）または式（１６）のようにいくつか定義することができる。
Figure 2021184148
Figure 2021184148
ここでγは実数の定数である。γが正の値であれば引力的な相互作用、γが負の値であれば斥力的な相互作用とみなせる。このＧ（ｘ）を用いて提案確率をｇ（ｘl→ｘl［ｊ0］）と与える。ｘl［ｊ0］はｘlのｊ0番目のビットをフリップした状態を表す。提案確率は、具体的には式（１７）または式（１８）のように定義することができる。
Figure 2021184148
Figure 2021184148
ここでＺ（ｘl）は規格化定数であり、計算方法は後述する。
レプリカ間距離としてハミング距離の１次関数を用いた場合、レプリカ間距離は式（１９）で定義できる。
Figure 2021184148
この場合、ΔＧ＝Ｇ（ｘl［ｊ0］）−Ｇ（ｘl）とｇ（ｘl→ｘl［ｊ0］）は以下のように計算できる。
Figure 2021184148
Figure 2021184148
このようにして、レプリカ間の相互作用を反映させた提案確率を計算することができる。次に、受け入れ確率の定義について説明する。
一般系の受け入れ確率ａ（ｘl→ｘl［ｊ0］）は、メトロポリス基準を採用することにすると、以下のように定義できる。
Figure 2021184148
すると遷移確率はＷ（ｘl→ｘl［ｊ0］）＝ｇ（ｘl→ｘl［ｊ0］）×ａ（ｘi→ｘi［ｊ0］）となる。従って、これらの計算に用いられる量はΔＥ、ΔＧ、規格化定数Ｚの３つである。
ここで、レプリカ間距離としてハミング距離の一次関数を用いた場合を例にとり、規格化定数Ｚの計算方法について説明する。規格化定数Ｚ（ｘl）は、提案候補がレプリカｌにおいて１つのビットをフリップしただけの状態なので、その総和として以下の式（２３）で計算される。
Figure 2021184148
式（２３）のまま規格化定数を計算しようとするとすべてのスピンの数だけ指数関数の和を計算することとなり、計算量が膨大となる。そこでイジングマシン３００は、１−ｂｉｔフリップであることに基づいて、差分計算（累積計算）を行うことで、計算量を抑止する。レプリカｌのｊ番目のビットだけをフリップした場合の規格化定数の差分は以下の式（２４）のようになる。
Figure 2021184148
イジングマシン３００は、式（２４）の右辺を計算することで求めた規格化定数の差分を、ビットフリップ前の規格化定数に加算することで、ビットフリップ後の規格化定数を求めることができる。なおイジングマシン３００は、ビットフリップが受け入れられた場合、そのときの規格化定数をレジスタまたはメモリに保存し、次回のビットフリップにおける規格化定数の算出に利用する。
次に、レプリカ間距離としてハミング距離の一次関数を用いた場合を例にとり、ΔＧの計算方法について説明する。ΔＧの計算は、一般にはレプリカ間の距離（あるいは距離の増加関数）の差分計算になる。単純な差分計算では、遷移前後でのレプリカ間のハミング距離を記憶しておかなければならない。しかし、具体的に距離（あるいは距離の増加関数）の形が分かれば、差分計算を行うことで、式（２５）、式（２６）に示すように、現在の状態だけに依存した量に書き換えることができる。
Figure 2021184148
Figure 2021184148
式（２６）において、ｘj0（ｌ）（ｘはチルダ付き）は、ｘj0（ｌ）（ｘはチルダ付き）＝（ｘj l，ｘj l，・・・，ｘj l）というビット列である。またｘj0（ｘはチルダ付き）は、ｘj0（ｘはチルダ付き）＝（ｘj l，ｘj 2，・・・，ｘj M）というビット列のベクトルである。
式（２６）を用いれば、レプリカ間距離をハミング距離の一次関数とした場合には、新たに導入したビット列のベクトル間のハミング距離だけで、ΔＧを記述することができる。従って、そのハミング距離だけを更新すればよい。
なお、これまでの説明では、１−ｂｉｔフリップの場合を想定しているが、１回の状態遷移で複数のビットをフリップする場合もある。例えばＯｎｅ−Ｈｏｔ制約がある問題を求解する場合である。
Ｏｎｅ−Ｈｏｔ制約は「ある変数の組の中で、値が１の変数は１個だけ」という制約である。この制約は、二次割り当て問題（ＱＡＰ：Quadratic Assignment Problem）、運搬経路問題（ＶＲＰ：Vehicle Routing Problem）などの様々な問題に適用される。
図７は、Ｏｎｅ−Ｈｏｔ制約下での１−ｂｉｔフリップの例を示す図である。図７の例では、イジングモデルの状態変数を示すビットが、４ビットずつのグループに分けられている。Ｏｎｅ−Ｈｏｔ制約では、同一グループに属するビットのうち１ビットのみが「１」であることを許される。このようなＯｎｅ−Ｈｏｔ制約の元で１−ｂｉｔフリップを行うと、１回の状態遷移では１ビットのみが反転し、制約違反の状態となる。１−ｂｉｔフリップをもう１回行うと、制約条件を満足することができる。
このようにＯｎｅ−Ｈｏｔ制約を有する問題を求解する場合、１−ｂｉｔフリップでは効率が悪い。そこで、イジングマシン３００は、１回の状態遷移で複数のビットをフリップさせることができる。
Ｏｎｅ−Ｈｏｔ制約には、１−Ｗａｙ １−Ｈｏｔ（１Ｗ１Ｈ）と２−Ｗａｙ １−Ｈｏｔ（２Ｗ１Ｈ）とがある。１Ｗ１Ｈは、１つの手段でビットをグルーピングしたときの各グループ内で、値が「１」のビットは１つのみという制約である。図７に示した例は１Ｗ１Ｈであり、１回の状態遷移で２つのビットをフリップさせることで、制約を満たしたままの状態遷移が可能となる。
２Ｗ１Ｈでは、２つの手段でビットがグルーピングされる。この場合、各ビットは、生成手段が異なる２つのグループに属する。そして２Ｗ１Ｈでも、各グループ内で値が「１」のビットは１つのみという制約がある。
図８は、２Ｗ１Ｈの制約を説明する図である。図８には、Ｎ個のビットをｎ×ｎ（ｎは１以上の整数）の正方行列の要素としている。Ｎ＝ｎ2である。２Ｗ１Ｈでは、各行および各列のビットの値の和は「１」であるという制約がある。すなわち同じ行のビットのうちの１つのビットだけ値が「１」であり、同じ列のビットのうちの１つのビットだけ値が１である場合に制約が満たされる。２Ｗ１Ｈの制約がある場合、１回の状態遷移で４つのビットをフリップさせることで、制約を満たしたままの状態遷移が可能となる。
ここでｍ＝１，２，・・・，Ｎとしたとき、１−ｂｉｔフリップ、１Ｗ１Ｈにおける２−ｂｉｔフリップ、２Ｗ１Ｈにおける４−ｂｉｔフリップそれぞれの、状態遷移、エネルギーの変化値ΔＥ、局所場の更新量Δｈは、以下のように表される。
＜１−ｂｉｔフリップ＞
・状態遷移 ：ｘi→ｘi＋Δｘi
・エネルギーの変化値：ΔＥi＝−ｈi・Δｘi
・局所場更新量 ：Δｈm＝Ｗm,i・Δｘi
＜１Ｗ１Ｈ（２−ｂｉｔフリップ）＞
・状態遷移：ｘi：１→０，ｘj：０→１
・エネルギーの変化値：ΔＥj＝ｈi−ｈj
・局所場更新量：Δｈm＝−Ｗm,i＋Ｗm,j
＜２Ｗ１Ｈ（４−ｂｉｔフリップ）＞
・状態遷移：ｘi：１→０，ｘj：０→１，ｘk：０→１，ｘl：１→０
・エネルギーの変化値：ΔＥj＝（ｈi＋ｈl）−（ｈj＋ｈk）−（Ｗil＋Ｗjk）
・局所場更新量：Δｈm＝Ｗmj＋Ｗmk−（Ｗmi＋Ｗml）
いずれの制約を適用するのかは、例えば、ユーザが問題の求解を指示する際に、ユーザによって指定される。イジングマシン３００は、指定された制約に応じたΔＥを計算し、レプリカ間の距離に応じた遷移確率で１または複数のビットを反転させる。
次に、レプリカ間距離を考慮したイジングマシン３００による解探索機能について説明する。
図９は、イジングマシンの解探索機能の一例を示す図である。イジングマシン３００は、データ受け取り部３４０、解探索エンジン３５０、および解出力部３６０を有している。データ受け取り部３４０と解出力部３６０とは、図４に示した制御回路３２０によって実現される機能である。解探索エンジン３５０は、図４に示した制御回路３２０が、ニューロン回路３１１，３１２，・・・，３１ｎとメモリ３３０とを制御することで実現する機能である。
データ受け取り部３４０は、制御装置２００から探索対象の問題の求解に用いる情報を受け取る。例えばデータ受け取り部３４０は、温度、レプリカ数、レプリカ相互作用の大きさ、iteration数（状態遷移の繰り返し回数）、初期状態などのパラメータを取得する。またデータ受け取り部３４０は、求解対象の問題を表すイジングモデルの重み係数を要素とする重み行列（２次式の係数）、バイアス行列（１次式の係数）、定数項、１−Ｈｏｔ制約のグループ情報などのデータを取得する。データ受け取り部３４０は、受け取った情報を解探索エンジン３５０に送信する。
解探索エンジン３５０は、複数のレプリカを用いて、エネルギーが最小となる解を探索する。そのために解探索エンジン３５０は、レプリカ保存部３５１と複数のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎとを有する。レプリカ保存部３５１は、例えば図４に示したメモリ３３０を利用して実現される。複数のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、イジングモデルに含まれるビットごとのニューロン回路を利用して実現される。
レプリカ保存部３５１は、レプリカの状態を記憶する。例えばレプリカが順番に更新されていくが、レプリカ間相互作用の計算には、更新前のレプリカの状態が使用される。そこでレプリカ保存部３５１が、更新前のレプリカの状態を記憶する。レプリカの状態は、状態変数に対応するビットの値、および温度パラメータなどのパラメータの値で表される。
各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、それぞれがレプリカによる解探索を行う。例えば各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、レプリカ保存部３５１を介して互いのレプリカの状態を示す情報をとやり取りしながらレプリカ間相互作用を計算し、解の探索を行う。
図１０は、解探索エンジンにおける処理の一例を示す図である。例えばレプリカ解探索部３５２ａは、重み係数（Ｗij）を記憶している。レプリカ解探索部３５２ａは、重み係数（Ｗij）と現在の各ビットの値（ｘ1 1,ｘ2 1，・・・，ｘN 1）とを用いて、式（４）に基づいて局所場（ｈ1，ｈ1，・・・，ｈN）を計算する。次にレプリカ解探索部３５２ａは、式（２６）に基づいて、各ビットがフリップした場合のレプリカ間の相互作用のエネルギーの差分（ΔＧ1，ΔＧ2，・・・，ΔＧN）を計算する。この際、レプリカ解探索部３５２ａは、レプリカ保存部３５１から他のレプリカの状態を示す情報（各ビットの値）を取得し、他のレプリカとの距離を計算し、その計算結果を用いてレプリカ間の相互作用のエネルギーの差分を計算する。
さらにレプリカ解探索部３５２ａは、局所場（ｈ1，ｈ2，・・・，ｈN）の値を用いて、エネルギーの変化値（Ｅ1，Ｅ2，・・・，ＥN）を計算する。なおエネルギーの変化値の計算式は、１−ｂｉｔフリップなのか１Ｗ１Ｈなのか２Ｗ１Ｈなのかによって異なる。例えば１−ｂｉｔフリップであれば、エネルギーの変化値は「ΔＥi＝−ｈi・Δｘi」である。１Ｗ１Ｈ（２−ｂｉｔフリップ）であれば、エネルギーの変化値は「ΔＥj＝ｈi−ｈj」である。２Ｗ１Ｈ（４−ｂｉｔフリップ）であれば、エネルギーの変化値はΔＥj＝（ｈi＋ｈl）−（ｈj＋ｈk）−（Ｗil＋Ｗjk）である。
レプリカ解探索部３５２ａは、エネルギーの変化値ΔＥから正のオフセット値Ｅoffを減算する。オフセット値Ｅoffには、フリップするビットが選択できなかった場合に、所定の値が加算される。オフセット値Ｅoffの増加は、フリップするビットが選択されるまで繰り返される。このように、オフセット値Ｅoffが増加することで、レプリカのエネルギーが極小値に留まる時間が短縮される。なお、オフセット値Ｅoffの初期値は、例えば「０」とする。
レプリカ解探索部３５２ａは、各ビットをフリップさせた場合のエネルギーの変化値ΔＥ（オフセット値Ｅoffが「０」以外の場合にはオフセット値Ｅoffを減算後の値）に基づいて、フリップするビット（更新ビット）を選択する。更新ビットの選択方法には、様々な方法がある（図１４〜図１７参照）。更新ビットの選択方法によっては、更新ビットの選択において、いずれのビットの更新の受け入れも棄却され、更新ビットが選択できないことがあり得る。レプリカ解探索部３５２ａは、例えば更新ビットが選択できなかった場合、オフセット値Ｅoffの値を増加させ、再度、更新ビットの選択を行う。
レプリカ解探索部３５２ａは、更新ビットが選択できた場合、更新ビットの値をフリップし、更新後のレプリカの状態「ｘ1 1，ｘ2 1，・・・，ｘN 1」を生成する。
レプリカ解探索部３５２ａ以外のレプリカ解探索部３５２ｂ，・・・，２５２ｎも、レプリカ解探索部３５２ａと同様に、更新後のレプリカの状態を生成する。
各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎが生成したレプリカの状態「ｘ1 1，ｘ2 1，・・・，ｘN 1」、「ｘ1 2，ｘ2 2，・・・，ｘN 2」、・・・、「ｘ1 N，ｘ2 N，・・・，ｘN N」は、レプリカ保存部３５１で保持される。各レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎは、レプリカ保存部３５１を参照することで、次回の状態更新時に、レプリカ間の相互作用のエネルギーの差分を算出することができる。
以下、解探索エンジン３５０による解探索の手順について詳細に説明する。
図１１は、解探索処理の手順の一例を示すフローチャートである。以下、図１１に示す処理をステップ番号に沿って説明する。
［ステップＳ１０１］解探索エンジン３５０は、複数のレプリカの初期状態（各ビットの値、温度パラメータの値など）を、そのレプリカの割り当て先のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎに設定する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカの初期状態に基づいて、初期エネルギー、初期のレプリカ間距離、初期規格化定数などを計算する。
［ステップＳ１０２］解探索エンジン３５０は、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎにレプリカごとの解探索を実行させる。レプリカごとの解探索処理の詳細は後述する（図１２参照）。
［ステップＳ１０３］解探索エンジン３５０は、解探索の終了条件を満たしたか否かを判断する。例えば解探索エンジン３５０は、ステップＳ１０２の処理の繰り返し回数が所定回数に達した場合に、終了条件を満たすと判断する。解探索エンジン３５０は、終了条件を満たした場合、処理をステップＳ１０８に進める。また解探索エンジン３５０は、終了条件が満たされていない場合、処理をステップＳ１０４に進める。
［ステップＳ１０４］解探索エンジン３５０は、複数のレプリカを温度パラメータの値で並べたときに隣接するレプリカの組を選択する。
［ステップＳ１０５］解探索エンジン３５０は、選択したレプリカの組の温度交換の実施の有無を決定する。例えば解探索エンジン３５０は、レプリカ間のエネルギーの差と各レプリカの温度パラメータの値とに基づいて、メトロポリスヘイスティング基準により交換確率を求める。そして解探索エンジン３５０は、交換確率が１であれば温度交換を実施すると判断する。また解探索エンジン３５０は、交換確率が１未満であれば、例えば０から１までの間の乱数を生成し、乱数の値が交換確率以下であれば、温度交換を実施すると判断する。
［ステップＳ１０６］解探索エンジン３５０は、温度交換を実施すると決定した場合、選択したレプリカの組それぞれの温度パラメータの値を交換する。
［ステップＳ１０７］解探索エンジン３５０は、隣接するレプリカのすべての組を選択したか否かを判断する。解探索エンジン３５０は、未選択の組がある場合、処理をステップＳ１０４に進める。また解探索エンジン３５０は、すべての組が選択済みの場合、処理をステップＳ１０２に進める。
［ステップＳ１０８］解探索エンジン３５０は、エネルギーが最小となるレプリカの状態を、解として出力する。
このようにして、レプリカ交換を行いながら、複数のレプリカを用いた効率的な解探索が行われる。
次にレプリカごとの解探索処理について詳細に説明する。
図１２は、レプリカごとの解探索処理の手順の一例を示すフローチャートである。以下、図１２に示す処理をステップ番号に沿って説明する。
［ステップＳ１２１］解探索エンジン３５０内のレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカについて、レプリカ間の相互作用のエネルギーの差分（ΔＧ1，ΔＧ2，・・・，ΔＧN）を計算する。レプリカ間の相互作用のエネルギーの差分の計算処理の詳細は後述する（図１３参照）。
［ステップＳ１２２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカについて、エネルギーの変化値（ΔＥ1，ΔＥ2，・・・，ΔＥN）を計算する。
［ステップＳ１２３］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、反復回数をインクリメントする。
［ステップＳ１２４］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、所定回数だけ反復したか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、所定回数だけ反復した場合、レプリカごとの解探索処理を終了する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、反復回数が所定回数に達していなければ、処理をステップＳ１２５に進める。
［ステップＳ１２５］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビット選択処理を行う。更新ビット選択処理の詳細は後述する（図１４〜図１７参照）。
［ステップＳ１２６］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択されたか否かを判断する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択されていない場合、処理をステップＳ１２５に進める。またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、更新ビットが選択された場合、処理をステップＳ１２７に進める。
［ステップＳ１２７］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカに関する情報を更新する。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、選択されたビットの状態をフリップさせ、各ビットの局所場ｈ、レプリカのエネルギーＥ、他のレプリカとのレプリカ間距離ｄ、規格化定数Ｚを更新する。その後、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは処理をステップＳ１２１に進める。
次に、レプリカ間の相互作用のエネルギーの差分（ΔＧ1，ΔＧ2，・・・，ΔＧN）の計算処理について詳細に説明する。
図１３は、レプリカ間の相互作用のエネルギーの差分の計算手順の一例を示すフローチャートである。以下、図１３に示す処理をステップ番号に沿って説明する。
［ステップＳ１４１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカと、そのレプリカ以外のレプリカそれぞれとのハミング距離を計算する。
［ステップＳ１４２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカのビットごとに、該当ビットをフリップした場合における遷移前後でのレプリカ間の相互作用のエネルギーの差分（ΔＧ1，ΔＧ2，・・・，ΔＧN）を計算する。例えば１番目のビットをフリップした場合のレプリカ間の相互作用のエネルギーの差分がΔＧ1である。
［ステップＳ１４３］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、割り当てられたレプリカの規格化定数Ｚを計算する。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカ間距離がハミング距離の一次式の場合は、状態遷移前後での規格化定数の差分を計算してもよい。差分を計算した場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、状態遷移ごとの規格化定数の差分を積算することで、最新の規格化定数を得ることができる。
次に、更新ビットの選択方法について説明する。更新ビットの選択方法としては、例えば以下の３つの方法が考えられる。
第１の更新ビット選択方法は、Original Boltzmannの方法である。第２の更新ビット選択方法は、エネルギーの並列計算を行い、エネルギーが下がる方向を先に参照することで効率的にビット更新を行う方法である。第３の更新ビット選択方法は、１イテレーションで常にビットフリップが起きるようにしたRejection-freeの方法である。
図１４は、第１の更新ビット選択方法による更新ビット選択処理の手順の一例を示すフローチャートである。以下、図１４に示す処理をステップ番号に沿って説明する。
［ステップＳ２０１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、レプリカ間距離を考慮に入れた提案確率ｇ（ｘl→ｘl［ｊ］）に従ってビットｊを選択する。
［ステップＳ２０２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、メトロポリス基準の受け入れ確率ａ（ｘl→ｘl［ｊ］）に従って、選ばれたビットをフリップするか否かを判定する。
第１の更新ビット選択方法は単純な方法であり、計算が容易であるが、選択したビットのフリップの提案が棄却されることもある。提案が棄却された場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、図１２のステップＳ１２６で「ＮＯ」と判断し、更新ビット選択処理を繰り返す。
第１の更新ビット選択方法は、提案分布に偏りがある影響を受けて受け入れ確率が小さくなってしまい、棄却ばかりが起こってしまう可能性がある。そこで更新ビットの提案が棄却された場合、レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、オフセット値Ｅoffの値を増加させることで、次回の更新ビットにおいて更新ビットが選択される確率を高めることができる。例えばレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、エネルギーが下がる方向が無くなる（エネルギー差がどのビット更新に対しても正になる）ときにはオフセット値Ｅoffに所定の値を加算する。
またレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、エネルギーの並列計算を行い、エネルギーが下がる方向を先に参照することで効率的にビット更新を行う第２の更新ビット選択方法を適用することもできる。
図１５は、第２の更新ビット選択方法の処理手順の一例を示すフローチャートである。以下、図１５に示す処理をステップ番号に沿って説明する。
［ステップＳ２１１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、すべてのビットに対して、メトロポリス基準の受け入れ確率ａ（ｘl→ｘl［ｊ］）に従って、該当ビットが選択された場合にフリップするか否かを判定する。レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットに対応付けて、判定結果を示すフラグを設定する。
［ステップＳ２１２］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットのフラグを参照し、ツリー状に接続されたセレクタを用いて、レプリカ間距離を考慮に入れた傾斜を与えて、更新ビットを選択する。
図１６は、更新ビット選択のためのツリー状に接続されたセレクタの一例を示す図である。制御回路３２０は、レプリカごとに、複数のビットそれぞれの状態遷移のエネルギーの変化値｛ΔＥｉ｝に応じて、上記の式（６）、式（７）の受け入れ確率でその状態遷移を許容するか否かを判断する。そして、制御回路３２０は、状態遷移を受け入れると判断したビットのうちの１つを、ツリー状に接続されたセレクタによって選択する。制御回路３２０は、選択したビットの番号と、遷移可否Ｆとを出力する。
このように、制御回路３２０は、複数のビットそれぞれに対して並列探索を行うことで、更新ビットが選択できる確率を高めることができる。
並列探索を行うため、制御回路３２０は、次の回路構成を有する。一例として、ビットの数を３２個として説明する。図１６の例ではいずれか１つのビットのみが更新ビットとして選択されるものとする。
制御回路３２０は、比較回路部５１〜５４とセレクタ部６０とを有する。
比較回路部５１〜５４は、複数の状態変数のそれぞれが遷移した場合のエネルギーの変化値｛ΔＥｉ｝を、ニューロン回路３１１，３１２，・・・，３１ｎから受け付ける。比較回路部５１〜５４は、｛ΔＥｉ｝に基づいて各状態遷移を受け入れるか否かを判定し、遷移可否｛ｆｉ｝を出力する。比較回路部５１〜５４それぞれは、８（＝３２／４）個の比較器を有する。比較回路部５１〜５４に含まれる全ての比較器の合計数は３２個となる。
具体的には、比較回路部５１は、比較器Ｃ０，Ｃ１，…，Ｃ７を有する。比較回路部５２は、比較器Ｃ８，Ｃ９，…，Ｃ１５を有する。比較回路部５３は、比較器Ｃ１６，Ｃ１７，…，Ｃ２３を有する。比較回路部５４は、比較器Ｃ２４，Ｃ２５，…，Ｃ３１を有する。比較器Ｃｉ（図１６の例ではｉは０以上３１以下の整数）は、ΔＥｉを受け付け、ΔＥｉに基づく判定に応じて受け入れ可否ｆｉを出力する。比較器Ｃｉによる判定では、エネルギーの変化値ΔＥｉと温度パラメータＴの値を用いて算出した受け入れ確率と、乱数値ｕ（０≦ｕ≦１）とが比較される。例えば比較器Ｃｉは、乱数値ｕが受け入れ確率以下であれば、ｉ番目のビットのフリップを受け入れると判定する。
比較回路部５１〜５４では、予め「Ｔ×ｌｏｇ（ｕ）」で表される値を計算することもできる。この値は、エネルギーが上がる状態遷移を確率的に引き起こす値であり、熱励起エネルギーまたは熱雑音と呼ぶこともできる。比較器Ｃｉは、ΔＥｉと熱励起エネルギーとを比較し、例えば熱励起エネルギーの方が大きければ、ｉ番目のビットのフリップを受け入れると判定する。
セレクタ部６０には、比較器Ｃｉの出力値が状態遷移の候補として入力される。そしてセレクタ部６０は、複数の状態遷移の候補の何れか１つを選択し、出力する。セレクタ部６０は、当該選択を行うためのｎ段（ｎは２以上の整数）のセレクタツリーを有する。図１６の例では、ｎ＝５である。
セレクタツリーの第１段目（１ｓｔ）は、部分セレクタ部６０ａ，６０ｂを有する。セレクタツリーの第２段目（２ｎｄ）は、部分セレクタ部６０ｃを有する。セレクタツリーの第３段目（３ｒｄ）は、部分セレクタ部６０ｄを有する。セレクタツリーの第４段目（４ｔｈ）は、部分セレクタ部６０ｅを有する。セレクタツリーの第５段目（５ｔｈ）は、部分セレクタ部６０ｆを有する。
部分セレクタ部６０ａ，６０ｂ，…，６０ｆのそれぞれは、例えば、２つの入力のうちの１つを選択用乱数により選択して出力する１または複数のセレクタを有する。セレクタ６１は、複数のセレクタのうちの１つであり、他のセレクタもセレクタ６１と同様の構成である。セレクタ６１に対する２つの入力はｉとｊの遷移番号を特定するための識別値Ｎｉ，Ｎｊと遷移可否情報ｆｉ，ｆｊと提案確率ｇ（ｘl→ｘl［ｉ］），ｇ（ｘl→ｘl［ｊ］）である。セレクタ６１の出力は遷移可否情報ｆｉ，ｆｊの論理和として得られる可否情報ｆｏと、ｉとｊのうち選択された方の遷移番号を特定するための識別値Ｎｏと、選択された方のビットの提案確率ｇ（ｘl→ｘl［ｏ］）である。
セレクタ６１は、遷移可否情報ｆｉ，ｆｊのいずれか一方が１（受け入れ可）、他方が０（受け入れ不可）の場合は受け入れ可の方のビットを必ず選択する。セレクタ６１は、遷移可否情報ｆｉ，ｆｊの両方０の場合はどのように選んでもよい。
セレクタ６１は、遷移可否情報ｆｉ，ｆｊの両方が１の場合には、候補選択用乱数を用いて、提案確率に応じた確率で一方を選択する。例えばセレクタ６１は、提案確率ｇ（ｘl→ｘl［ｉ］），ｇ（ｘl→ｘl［ｊ］）の比率に応じて、０から１の値域を、ｉとｊのビットに対応する２つの区間に分ける。そしてセレクタ６１は、候補選択用乱数を含む区間に対応するビットを選択する。そしてセレクタ６１は、選択結果により選ばれたビットの識別値Ｎｏを生成し出力する。
図１６の例ではセレクタ６１以外のセレクタが略記されている。図１６では黒い丸印で表された箇所が、１つのセレクタに相当する。部分セレクタ部６０ａ，６０ｂ，６０ｃのそれぞれは、８個のセレクタを有する。部分セレクタ部６０ｄは、４個のセレクタを有する。部分セレクタ部６０ｅは、２個のセレクタを有する。部分セレクタ部６０ｆは、１個のセレクタを有する。部分セレクタ部６０ａ〜６０ｆ内の各セレクタがセレクタ６１と同様の選択処理を行うことで、レプリカ間距離に応じた提案確率が高いビットほど選択される可能性を高くして、１つのビットが状態遷移の候補として出力される。
図１６で示されるように、制御回路３２０は、状態遷移の並列探索を行い、セレクタの２進木構造を用いてノックダウン方式（あるいはトーナメント方式とも呼ばれる）で、状態遷移の候補を１つに絞り込む。フリップによりエネルギーが減少するビットは、比較器により受け入れ可と判断されるため、フリップによりエネルギーが減少するビットが少なくとも１つ存在すれば、セレクタ部６０による１回の選択で更新ビットを選択できる。また局所解に達しており、いずれのビットをフリップしてもエネルギーが増加する場合であっても、乱数値ｕと温度パラメータＴの値に基づいて、いずれか１つのビットのフリップが受け入れられる可能性がある。いずれか１つのビットのフリップが受け入れられれば、セレクタ部６０による１回の選択で更新ビットを選択できる。しかも、セレクタによる選択時に、レプリカ間距離を反映させた提案確率を用いたことで、提案確率が高いビットほど更新ビットとして選択される可能性が高くなる。
なおレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、セレクタ部６０が出力した遷移可否情報が０（受け入れ不可）の場合、オフセット値を増加させて、更新ビット選択処理を繰り返す。これにより、更新ビットを早期に選択できる可能性が高くなる。
図１７は、第３の更新ビット選択方法の処理手順の一例を示すフローチャートである。第３の更新ビット選択方法は、以下の１ステップで更新ビットを選択できる。
［ステップＳ２３１］レプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、各ビットの遷移確率Ｗ（ｘl→ｘl［ｊ0］）＝ｇ（ｘl→ｘl［ｊ0］）×ａ（ｘl→ｘl［ｊ0］を用いて、以下の式（２７）に示すRejection-freeの遷移確率Ｗ（ｘl→ｘl［ｊ0］）（Ｗはチルダ付き）を計算する。
Figure 2021184148
そしてレプリカ解探索部３５２ａ，３５２ｂ，・・・，３５２ｎそれぞれは、Ｒｅｊｅｃｔｉｏｎ−ｆｒｅｅの遷移確率により、いずれか１つのビットを更新ビットとして選択する。このように各ビットの遷移確率を正規化し、受け入れ確率の合計が１となるようにすることで、１回の更新ビット選択処理で、常に更新ビットを選択することが可能となる。
以上説明したように、第２の実施の形態に係るイジングマシン３００は、レプリカ間の相互作用を提案確率に反映させ、複数のレプリカを用いた解探索を行っている。これにより、組み合わせ最適化問題をメトロポリスヘイスティングの方法に基づいて求解する際に、収束先の分布を保ったまま、それぞれのレプリカがバラバラに状態空間を探索することが期待され、求解性能が向上する。すなわち、最適解へ到達する可能性が高くなり、エネルギーの下がり方を速くすることができる。
図１８は、レプリカ間に斥力の相互作用を設定した場合のエネルギーランドスケープを示す図である。複数のレプリカ７１〜７３の間には、斥力の相互作用が設けられている。この場合、レプリカ７１〜７３が互いに反発し合うことで、広い探索空間を効率的に探索することができる。例えば図１８の例では、複数のレプリカ７１〜７３それぞれが、異なる探索範囲を探索している。このような複数の探索範囲が複数のレプリカ７１〜７３で分担して探索されることで、広範囲を効率的に探索することができる。
図１９は、レプリカ間に引力の相互作用を設定した場合のエネルギーランドスケープを示す図である。複数のレプリカ７４，７５の間には、引力の相互作用が設けられている。レプリカ７４，７５が互いに引きつけられることで局所解から脱出し易くなり、集団全体として大域解に到達できる可能性を高まる。
次に、効果確認を行った検証例について、図２０と図２１を参照して説明する。図２０は、第１の検証例を示す図である。図２１は、第２の検証例を示す図である。図２０、図２１に示す例は、二次割り当て問題（ＱＡＰ）という代表的な組み合わせ最適化問題のいくつかのインスタンスについて検証した結果である。提案分布に従った各ビットの提案確率の計算には前述の式（１７）を使用している。レプリカ間の相互作用のエネルギーとしては、前述の式（１９）に示したハミング距離の一次関数を使用している。更新ビットの選択方法としては、第３の更新ビット選択方法（Rejection-free）が用いられいている。
図２０、図２１の例では、１−ｂｉｔフリップ遷移かつレプリカ交換を使う解探索手法において、レプリカ間の相互作用の有無によるエネルギーの下がり方の違いを比較している。横軸が状態遷移の反復回数であり、縦軸がその時点で得られているエネルギーの最小値である。γ（図２０、図２１では“ｇａｍｍａ”と表記）を斥力相互作用のパラメータとしたときに、ｇａｍｍａ＝０とｇａｍｍａ＜０（つまり斥力相互作用の有無）との場合についてエネルギーの下がり方について比較している。
図２０の例では、斥力相互作用を導入した場合（ｇａｍｍａ−３）のほうが相互作用を導入しない場合（ｇａｍｍａ−０）よりもエネルギーの下がり方が速い。図２１の例においても、斥力相互作用を導入した場合（ｇａｍｍａ−１０）のほうが相互作用を導入しない場合（ｇａｍｍａ−０）よりもエネルギーの下がり方が速い。
このように、レプリカ間の相互作用を導入したことで、解探索性能が向上している。しかもレプリカ間の相互作用を提案確率に反映しており、目的関数に手を加えないため、適切な目的関数（例えばギブス分布を示す関数）を用いた解探索が可能となる。
なお、非特許文献１に示されたＣＭＣと呼ばれる方法は、実数を定義域とする目的関数にのみ適用可能な方法であり、２値の離散空間を定義域とする（バイナリ変数）イジングマシンの目的関数に直接は適応できない。またＣＭＣでは、距離が近いレプリカの数（密度）をカウントしているが、１−ｂｉｔフリップの場合にレプリカすべての状態をみたときに、その状況がフリップした前後で大きくは変わらない。そのため、あるビットのフリップ前後でのレプリカ数の密度の比はほぼ１に近くなってしまい、２値の離散空間を定義域とするとレプリカ相互作用の効果が薄くなってしまう。それに対して、第２の実施の形態に示した方法では、２値の離散空間を定義域とする組み合わせ最適化問題に適用でき、求解性能も向上する。
また、非特許文献２に示されたＲＥと呼ばれる方法では、レプリカ間の相互作用を目的関数に直接加える方法をとっているため、本来の目的関数の最適化を行っている保証はない。それに対して、第２の実施の形態に示した方法では、レプリカ間の相互作用を提案確率に反映しており、適切な目的関数を用いた解探索が可能となる。
〔その他の実施の形態〕
第２の実施の形態では、レプリカ間の温度交換を行っているが、レプリカ間の温度交換を実施せずに、複数のレプリカで個別に解探索を行うことも可能である。その場合であっても、レプリカ間の相互作用を考慮した解探索により、解探索能力が向上する。
また第２の実施の形態では、２値の離散空間を定義域とするイジングモデルを用いた求解を行っているが、実数を定義域とするモデルをレプリカとして求解する場合にも適用可能である。
さらに第２の実施の形態では、複数のニューロン回路３１１，３１２，・・・，３１ｎを有するイジングマシン３００で解探索を行っているが、同じ処理を図２に示したサーバ１００と同様のハードウェア構成のノイマン型コンピュータで実現することも可能である。その場合、イジングマシン３００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態と同様の解探索処理を実行する。イジングマシン３００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、イジングマシン３００に実行させるプログラムをストレージ装置に格納しておくことができる。イジングマシン３００のプロセッサは、ストレージ装置内のプログラムの少なくとも一部をメモリにロードし、プログラムを実行する。またイジングマシン３００に実行させるプログラムを、光ディスク、メモリ装置、メモリカードなどの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばイジングマシン３００のプロセッサからの制御により、ストレージ装置にインストールされた後、実行可能となる。またイジングマシン３００のプロセッサが、可搬型記録媒体から直接プログラムを読み出して実行することもできる。
以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。
１ 状態空間
２〜４ レプリカ
１０ 最適化装置
１１ 記憶部
１２ 処理部
Figure 2021184148
右辺の１項目は、Ｎ個の状態変数の全組み合わせについて、漏れと重複なく、２つの状態変数の値（０または１）と重み係数との積を積算したものである。ｘiはｉ番目の状態変数、ｘjはｊ番目の状態変数を表し、Ｗijは、ｘiとｘjとの結合の強さを示す重み係数である。右辺の２項目は、各状態変数のそれぞれについてのバイアス係数（ｂi）とｘ ｉ との積の総和を求めたものである。Ｗijが正の場合、ｘiとｘjとが同じ値となるような相互作用が働く。またＷijが負の場合、ｘiとｘjとが異なる値となるような相互作用が働く。なおＷij＝Ｗjiであり、Ｗii＝０である。

Claims (13)

  1. 複数のレプリカそれぞれの複数の状態変数の値を記憶する記憶部と、
    前記複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での前記レプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する処理部と、
    を有する最適化装置。
  2. 前記処理部は、さらに、
    前記第１状態変数の値を更新すると決定した場合、前記第１状態変数の値を更新した後の前記レプリカの前記複数の状態変数の値の組み合わせに基づいて、目的関数の値を特定すると共に、前記記憶部に記憶された前記レプリカの前記第１状態変数の値を更新し、
    前記複数のレプリカそれぞれの前記複数の状態変数の一の状態変数の値の更新を繰り返し、前記目的関数の値が所定の条件を満たしたときの前記複数の状態変数の値の組み合わせを出力する、
    請求項１記載の最適化装置。
  3. 前記処理部は、前記第１状態変数の値を更新すると前記レプリカと前記他のレプリカとの前記距離が遠ざかる場合に前記相互作用の強さを増加させ、前記相互作用の強さの増加量が大きいほど前記提案確率を大きくする、
    請求項１または２記載の最適化装置。
  4. 前記処理部は、前記第１状態変数の値を更新すると前記レプリカと前記他のレプリカとの前記距離が近づく場合に前記相互作用の強さを増加させ、前記相互作用の強さの増加量が大きいほど前記提案確率を大きくする、
    請求項１または２記載の最適化装置。
  5. 前記処理部は、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量をΔＧ、前記レプリカに設定された温度パラメータの値の逆数である逆温度をβとしたとき、ｅｘｐ（−βΔＧ）を所定の規格化定数で除算した値を前記提案確率とする、
    請求項１ないし４のいずれかに記載の最適化装置。
  6. 前記処理部は、前記第１状態変数の値を更新した場合における相互作用の強さの前記変化量をΔＧ、前記レプリカに設定された温度パラメータの値の逆数である逆温度をβとしたとき、１とｅｘｐ（−βΔＧ）とのうちの小さい方を、所定の規格化定数で除算した値を前記提案確率とする、
    請求項１ないし４のいずれかに記載の最適化装置。
  7. 前記処理部は、前記複数の状態変数それぞれを前記第１状態変数とした場合における前記複数の状態変数それぞれについてのｅｘｐ（−βΔＧ）の値の総和を、前記規格化定数とする、
    請求項５または６記載の最適化装置。
  8. 前記処理部は、前記レプリカが更新する前記第１状態変数の決定と、決定された前記第１状態変数の更新とを繰り返し、前記規格化定数の計算では、前回に更新された前記第１状態変数の決定の際に用いられた前記規格化定数の値と、前回の前記第１状態変数の更新の前と更新後との前記規格化定数の値の差分とに基づいて、今回に更新する前記第１状態変数の決定で用いる前記規格化定数の値を算出する、
    請求項５ないし７のいずれかに記載の最適化装置。
  9. 前記処理部は、前記他のレプリカそれぞれとの前記距離の合計に基づく値を、前記相互作用の強さとする、
    請求項１ないし８のいずれかに記載の最適化装置。
  10. 前記処理部は、前記他のレプリカそれぞれとの前記距離の平方根の合計に基づく値を、前記相互作用の強さとする、
    請求項１ないし８のいずれかに記載の最適化装置。
  11. 前記処理部は、前記複数の状態変数それぞれについて、前記受け入れ確率に基づいて、状態変数の更新が提案された場合に更新を受け入れるか否かを決定し、更新を受け入れると判定された状態変数の中から、前記提案確率が高い状態変数ほど選択される可能性を高くして、少なくとも１つの状態変数を更新対象に決定する、
    請求項１ないし１０のいずれかに記載の最適化装置。
  12. 最適化装置が、
    複数の状態変数を有する複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での前記レプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を特定し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する、
    を有する最適化方法。
  13. 最適化装置に、
    複数の状態変数を有する複数のレプリカそれぞれについて、該レプリカが有する前記複数の状態変数のうちの第１状態変数の値を更新した場合における、前記複数の状態変数の値の組み合わせが取り得る空間を示す状態空間内での前記レプリカと他のレプリカとの距離の変化に応じた相互作用の強さの変化量を計算し、前記第１状態変数の値を更新した場合における前記相互作用の強さの前記変化量に応じた提案確率と、目的の確率分布に応じた受け入れ確率とに基づいて、前記第１状態変数の値を更新するか否かを決定する、
    処理を実行させる最適化プログラム。
JP2020088882A 2020-05-21 2020-05-21 最適化装置、最適化方法、および最適化プログラム Withdrawn JP2021184148A (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2020088882A JP2021184148A (ja) 2020-05-21 2020-05-21 最適化装置、最適化方法、および最適化プログラム
EP21153978.8A EP3913506A1 (en) 2020-05-21 2021-01-28 Optimization device, optimization method, and optimization program
US17/168,557 US11790130B2 (en) 2020-05-21 2021-02-05 Optimization device, optimization method, and non-transitory computer-readable storage medium for storing optimization program
CN202110183741.8A CN113705851A (zh) 2020-05-21 2021-02-10 优化装置、优化方法和用于存储优化程序的存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2020088882A JP2021184148A (ja) 2020-05-21 2020-05-21 最適化装置、最適化方法、および最適化プログラム

Publications (1)

Publication Number Publication Date
JP2021184148A true JP2021184148A (ja) 2021-12-02

Family

ID=74346946

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2020088882A Withdrawn JP2021184148A (ja) 2020-05-21 2020-05-21 最適化装置、最適化方法、および最適化プログラム

Country Status (4)

Country Link
US (1) US11790130B2 (ja)
EP (1) EP3913506A1 (ja)
JP (1) JP2021184148A (ja)
CN (1) CN113705851A (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023228336A1 (ja) * 2022-05-25 2023-11-30 Tdk株式会社 計算モデル及び計算プログラム

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003223322A (ja) * 2002-01-30 2003-08-08 Mitsubishi Electric Corp 組合せ最適化問題の解析装置
JP5865457B1 (ja) * 2014-08-29 2016-02-17 株式会社日立製作所 情報処理システム及び管理装置
JP6468254B2 (ja) 2016-07-01 2019-02-13 富士通株式会社 情報処理装置、イジング装置及び情報処理装置の制御方法
JP6979331B2 (ja) 2017-10-30 2021-12-15 株式会社日立製作所 情報処理装置および情報処理方法
JP7007585B2 (ja) * 2018-03-16 2022-01-24 富士通株式会社 最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
JP7004906B2 (ja) 2018-05-01 2022-01-21 富士通株式会社 最適化装置及び最適化装置の制御方法
EP3852029A4 (en) * 2018-09-13 2022-06-15 Fujitsu Limited OPTIMIZATION DEVICE AND CONTROL METHOD OF AN OPTIMIZATION DEVICE
JP7206476B2 (ja) * 2018-09-14 2023-01-18 富士通株式会社 最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
JP7193707B2 (ja) * 2018-10-11 2022-12-21 富士通株式会社 最適化装置、最適化装置の制御方法、サンプリング装置およびサンプリング装置の制御方法
JP2019071119A (ja) 2019-01-11 2019-05-09 富士通株式会社 情報処理装置、イジング装置及び情報処理装置の制御方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023228336A1 (ja) * 2022-05-25 2023-11-30 Tdk株式会社 計算モデル及び計算プログラム

Also Published As

Publication number Publication date
CN113705851A (zh) 2021-11-26
EP3913506A1 (en) 2021-11-24
US11790130B2 (en) 2023-10-17
US20210365605A1 (en) 2021-11-25

Similar Documents

Publication Publication Date Title
Wang et al. Quantumnas: Noise-adaptive search for robust quantum circuits
Mirhoseini et al. A graph placement methodology for fast chip design
Lujan-Moreno et al. Design of experiments and response surface methodology to tune machine learning hyperparameters, with a random forest case-study
Devezer et al. Scientific discovery in a model-centric framework: Reproducibility, innovation, and epistemic diversity
Alazzam et al. A new optimization algorithm for combinatorial problems
Dong et al. Adaptive neural network-based approximation to accelerate eulerian fluid simulation
JP7488457B2 (ja) 最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
Nadimi-Shahraki et al. MMKE: Multi-trial vector-based monkey king evolution algorithm and its applications for engineering optimization problems
Pandey et al. Evaluation of genetic algorithm’s selection methods
JP6925546B1 (ja) 演算システム、情報処理装置、および最適解探索処理方法
JP2022090249A (ja) 最適化装置、最適化方法、及び最適化プログラム
Gao et al. A novel hybrid PSO based on levy flight and wavelet mutation for global optimization
JP2021005282A (ja) 最適化装置、最適化装置の制御方法及び最適化装置の制御プログラム
Khan Particle swarm optimisation based feature selection for software effort prediction using supervised machine learning and ensemble methods: A comparative study
JP2021184148A (ja) 最適化装置、最適化方法、および最適化プログラム
JP2022150078A (ja) 情報処理プログラム、情報処理装置、及び情報処理方法
Steffen Particle swarm optimization with a simplex strategy to avoid getting stuck on local optimum
JP7215966B2 (ja) ハイパーパラメータ管理装置、ハイパーパラメータ管理方法及びハイパーパラメータ管理プログラム製品
Kang et al. TMHSCA: a novel hybrid two-stage mutation with a sine cosine algorithm for discounted {0-1} knapsack problems
CN111582313A (zh) 样本数据生成方法、装置及电子设备
JPWO2019208564A1 (ja) ニューラルネットワーク学習装置、ニューラルネットワーク学習方法、プログラム
Briffoteaux Parallel surrogate-based algorithms for solving expensive optimization problems
Zhu et al. A hybrid model for nonlinear regression with missing data using quasilinear kernel
JP2022052222A (ja) 最適化装置、最適化方法、および最適化プログラム
CN116010754A (zh) 存储程序的计算机可读记录介质、数据处理方法和设备

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20200923

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20230209

A761 Written withdrawal of application

Free format text: JAPANESE INTERMEDIATE CODE: A761

Effective date: 20240129