(57)【要約】【課題】より精度のよい最適解を求める。【解決手段】実施形態の最適化装置は、変換部と、２進展開部と、イジング処理部と、逆２進展開部とを有する。変換部は、入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換する。２進展開部は、第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする。イジング処理部は、第２のコスト関数の係数に関するイジング形式の値について、基底状態を示す値の探索を行う。逆２進展開部は、探索した基底状態を示す値について逆２進展開を行い、第１のコスト関数を最小とする値を求める。【選択図】図１
2020067811.tif 000002
【特許請求の範囲】
【請求項１】
入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換する変換部と、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする２進展開部と、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行うイジング処理部と、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める逆２進展開部と、
を有することを特徴とする最適化装置。
【請求項２】
前記変換部は、前記第１のコスト関数に含まれる不連続関数について、ルジャンドル変換を行い、当該ルジャンドル変換後の式をもとにＷｏｌｆｅの双対定理による双対問題を求めて得られた関数を、前記不連続関数に代えて前記第１のコスト関数に盛り込むことで、前記第２のコスト関数に変換する、
ことを特徴とする請求項１に記載の最適化装置。
【請求項３】
入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換し、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とし、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行い、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める、
処理をコンピュータが実行することを特徴とする最適化方法。
【請求項４】
入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換し、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とし、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行い、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める、
処理をコンピュータに実行させることを特徴とする最適化プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、最適化装置、最適化方法および最適化プログラムに関する。
【背景技術】
【０００２】
従来、ノイマン型コンピュータが不得意とする多変数の最適化問題を解く方法として、イジング型のエネルギー関数（コスト関数、目的関数とも呼ばれる）を用いたイジングマシン（ボルツマンマシンと呼ばれる場合もある）がある。イジングマシンでは、計算対象の問題を、磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えて計算する。例えば、イジングマシンは、シミュレーテッド・アニーリング（ＳＡとも呼ばれる）により、上記のようなエネルギー関数の最小値が得られる各ビットの状態の組み合わせ（基底状態）を、解として求める。
【０００３】
計算対象の問題のイジングモデルへの置き換えについては、問題とするコスト関数に不連続関数の項が存在するとイジングモデルへの変換ができないことから、ｑ−ｌｏｓｓ関数で近似してイジングモデルに置き換える従来技術が知られている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開平１０−１４９４４５号公報
【特許文献２】特表２００２−５２２１２８号公報
【非特許文献】
【０００５】
【非特許文献１】V.S.  Denchev,  N.  Ding,  S.V.N.  Vishwanathan,  and  H.  Neven,  “Robust  classication  with  adiabatic  quantum  optimization,“  in  Proc.  ICML'12,  pp.  1003-1010  (2012).
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上記の従来技術では、不連続関数をｑ−ｌｏｓｓ関数で近似した別関数に置き換えることから、イジングマシンにより求めた解と、本来求めたい値との間に差が生じ、高精度な解が得られない場合があるという問題がある。
【０００７】
図１０は、従来の最適化装置の機能構成例を示すブロック図である。図１０に示すように、最適化装置２００では、最適化問題に関するコスト関数Ｆ（Ｗ）の入力を受け付け、ｑ−ｌｏｓｓ関数を用いてイジングモデルの式（Ｆ’（Ｗ））への変換を行う（２０１）。最適化問題に関するコスト関数Ｆ（Ｗ）としては、一例として、医療分野における放射線照射計画の策定に関するものがある。
【０００８】
図１１は、放射線照射計画の策定例を説明する説明図である。図１１に示すように、ＩＭＲＴ（強度変調放射線治療）においては、マルチリーフコリメータ３０３の形状を変化させ、マルチリーフコリメータ３０３により遮られていない部分の形状に合わせて細かい解像度で強度を変化させた放射線を腫瘍３０１および正常組織３０２などの部位ごとに照射する。この際の放射線照射計画は、正常組織３０２への照射を低減しつつ、腫瘍３０１への照射を行うような最適化問題として捉えることができる。
【０００９】
この放射線照射計画に関するコスト関数（Ｆ（Ｗ））は、次のとおりである。ここで、図１１における照射対象（腫瘍３０１、正常組織３０２）は、２次元で説明しており、以下では１次元に展開しているので、照射対象が３次元以上あっても以降は同様であるものとする。
【００１０】
まず、照射対象を離散化（図示例における四角での区切り）し、一列に並べたピクセルとする。ここで、ピクセルの個数はＩ個とする。また、放射線のビームはＮ本あるものとし、その強さをベクトルＷ＝（ｗ１，Ｗ２，…，ｗＮ）で表す。
【００１１】
ｉ番目のピクセルにおける放射線の照射量をｄｉ（Ｗ，Ｍｉ）とすると、ｄｉ（Ｗ，Ｍｉ）は次の式（１）のとおりとなる。
【００１２】
【数１】
2020067811.tif 000003
【００１３】
ただし、Ｍｉｎは、ｎ番目のビームがｉ番目のピクセルに与える影響であり、事前に計算可能である。また、Ｍｉ＝（Ｍｉ１，Ｍｉ２，…，ＭｉＮ）とする。
【００１４】
ピクセルごとに目的の照射量（ｐｉ）が与えられた場合、これに関するコスト関数Ｆ（Ｗ）は、目標とする照射量をＭ＝（Ｍ１，Ｍ２，…，Ｍｉ）として次の式（２）のとおりとなる。なお、式（２）において二乗しているのは、イジングモデルに変換するためと、ｐｉとｄｉ（Ｗ，Ｍｉ）が最も近い値となったときに最小となるＷを求める際に差分の符号を無視したいためである。
【００１５】
【数２】
2020067811.tif 000004
【００１６】
しかし、臓器によっては照射量に限界がありこれをＣｉと置くと、コスト関数Ｆ（Ｗ）は次の式（３）のようになる。
【００１７】
【数３】
2020067811.tif 000005
【００１８】
この式（３）に示すコスト関数Ｆ（Ｗ）を最小にするビームの強さのベクトルＷが、最適化問題として求める答えとなる。ここで、ｐｉ，Ｍｉ，βｉ、Ｃｉは、ユーザーにより設定される定数である。また、式（３）に示すコスト関数Ｆ（Ｗ）における第１項は、連続関数となる項（Ｃ（Ｗ））である。第２項は、制約項（Ｆ１（Ｗ））であり、βｉはピクセルごとに制約項の強さをあらわす。ｄｉ（Ｗ）がＣｉを超えて大きくなるほどＦ（Ｗ）も大きくなる。また、第２項が不連続なステップ関数となるため、Ｆ（Ｗ）はこのままではイジング形式で表現できない。
【００１９】
そこで、コスト関数Ｆ（Ｗ）をイジングモデルの式（Ｆ’（Ｗ））へ変換する処理（２０１）では、Ｆ（Ｗ）に含まれるＦ１（Ｗ）を二乗してｑ−ｌｏｓｓ関数に置き換え、２体イジング形式で表現したＦ’（Ｗ）を出力する。
【００２０】
図１２は、イジングモデルの式への変換の従来例を示すフローチャートである。図１２に示すように、処理が開始されると、最適化装置２００は、式（３）に示すコスト関数Ｆ（Ｗ）の入力を受け付ける（Ｓ２０１）。
【００２１】
このコスト関数Ｆ（Ｗ）には、不連続関数となるＦ１（Ｗ）が含まれる。また、ピクセル数（Ｉ）、ビームの本数（Ｎ）、目標照射量（ｐｉ）、ビームが与える影響度（Ｍｉ）、制約の強さ（βｉ）、照射限界量（Ｃｉ）は、コスト関数Ｆ（Ｗ）に定数として含まれる。この定数はユーザーが指定する。
【００２２】
次いで、最適化装置２００は、ｑ−ｌｏｓｓ関数で処理するために、コスト関数Ｆ（Ｗ）に含まれる不連続関数（Ｆ１（Ｗ））を二乗する（Ｓ２０２）。
【００２３】
次いで、最適化装置２００は、不連続関数（Ｆ１（Ｗ））を二乗したコスト関数Ｆ（Ｗ）をｑ−ｌｏｓｓ関数を使ってイジング形式に変換する（Ｓ２０３）。次いで、最適化装置２００は、変換したＦ’（Ｗ）を出力する（Ｓ２０４）。
【００２４】
ｑ−ｌｏｓｓ関数は、まず、パラメータ（ｑ∈（−∞，０]）を使って次の式（４）のように書ける。
【００２５】
【数４】
2020067811.tif 000006
【００２６】
ここで、変数（ｔ∈Ｒ）を加えてイジング形式に書き直すと、次の式（５）のとおりとなる。
【００２７】
【数５】
2020067811.tif 000007
【００２８】
ここで、ｑ＜＜０として、式（４）はＬｑ（ｍ）＝（ｍａｘ［０，１−ｍ］）２となるため、Ｆに含まれるＦ１はＬｑを用いて次の式（６）のとおりとなる。
【００２９】
【数６】
2020067811.tif 000008
【００３０】
また、イジング形式では、ｑ＜＜０と、最小条件であるｔｉ≧１の条件をつけると、次の式（７）となる。
【００３１】
【数７】
2020067811.tif 000009
【００３２】
この式（７）において、ｔを定数１にして（定数項１と相殺する）次の式（８）のＦ’（Ｗ）を得る。このＦ’（Ｗ）では、制約項は二乗の形であらわされている。最適化装置２００は、得られたＦ’（Ｗ）を出力する。
【００３３】
【数８】
2020067811.tif 000010
【００３４】
図１０に戻り、最適化装置２００は、変換したＦ’（Ｗ）が連続値のままであることから、イジングマシンで処理するために、Ｗを０，１に２進展開する（２０２）。具体的には、最適化装置２００は、２進展開の精度（ｄｗ）、値の範囲（α、β）の入力を受け付け、Ｗの２進展開を行う。次いで、最適化装置２００は、２進展開を行ったＷ’に関する結合整数（ｑ）およびＷ’を元の実数に戻すための係数（Ｚ）を出力する。
【００３５】
図１３は、２進展開処理の一例を示すフローチャートである。図１３に示すように、処理が開始されると、最適化装置２００は、変換したイジング形式の式（Ｆ’（Ｗ））と、Ｗの精度（ビット深度）ｄｗと、Ｗの値の範囲αｗ〜βｗとの入力を受け付ける（Ｓ３０１）。ここで、ｄｗ、αｗ〜βｗは、ユーザーが与える数値である。
【００３６】
次いで、最適化装置２００は、２進展開した結果得られる各ビットの重み（δ）を計算し、δｗ（ｋ）＝２ｋ−１／（２ｋ−１−１）とする（Ｓ３０２）。ただし、ｋ＝１〜ｄｗである。
【００３７】
例えば、数値Ｗを２進数にしたものに対してビット精度（ｄｄ）で２進展開する際に、ＭＳＢを左側としてｄｄビットまで２進数にした配列ｗｄをｗｎｄ＝［ｗｄ１、ｗｄ２，…，ｗｄｄｄ］とする。これに対応する重みδｗは、展開したビットと重みを乗算して次の式（９）となる。
【００３８】
【数９】
2020067811.tif 000011
【００３９】
次いで、最適化装置２００は、Ｗｎの区間をαｗ〜βｗとして、２進展開で同時に範囲を抑えると、Ｗｎは次の式（１０）のとおりとなる。
【００４０】
【数１０】
2020067811.tif 000012
【００４１】
すなわち、最適化装置２００は、ｎ＝１〜Ｎとして、Ｗｎ＝αｗΣｋ＝１ｄｗＷｎ，ｋδｗ（ｋ）＋βｗとする（Ｓ３０３）。なお、Ｗｎ，ｋは、Ｗｎを２進数にしたもので、０，１の数値であり、イジングスピンに対応する。
【００４２】
次いで、最適化装置２００は、上記式をＦ（ｗ，ｔ）に代入してえた式（Ｚ）のＷｎの２次の項と１次の項の係数を結合係数として出力する（Ｓ３０４）。このとき、最適化装置２００は、Ｚ＝［αｗ，βｗ，δｗ］も出力する。
【００４３】
具体的には、Ｆ’のすべての項はｉについて総和しているため最後に足し合わせるものとし、まず総和の内部の項を考える。ここで、Ａ＝Σｎ＝１ＮＷｎＭｎとするとＦ’の内部は次の式（１１）のとおりとなる。
【００４４】
【数１１】
2020067811.tif 000013
【００４５】
ここで、Ｓｒ＝［Ｗ１，Ｗ２，…，ＷＮ］、Ｍ＝［Ｍ１，Ｍ２，…，ＭＮ］とする。
【００４６】
式（１１）において、（１＋β）Ａ２は、２次の項であり、以下の結合係数（Ｊｒ）に対応する。また、−２（Ｐ＋βＣ）Ａは、１次の項であり、以下の結合係数ｈｒに対応する。また、ｐ２＋βＣ２は、定数の項であり、Ｗと無関係なので以下では無視する。
【００４７】
２進展開する前のイジングモデルの式（Ｓはイジングスピンで実数）と結合係数Ｊｒ、ｈｒを次の式（１２）のとおりである。
【００４８】
【数１２】
2020067811.tif 000014
【００４９】
よって、結合係数Ｊｒ、ｈｒは、次の式（１３）、（１４）のとおりとなる。
【００５０】
【数１３】
2020067811.tif 000015
【００５１】
【数１４】
2020067811.tif 000016
【００５２】
Ｓｒは、実数なので２進展開する。２進展開はＳｒの各要素を行方向に２進展開した要素を並べたものである。ここで、２進展開したＳｒをＳ＝［Ｗ１，１，Ｗ１，２，…，ＷＮ，ｄ］とする。重みδｗをＪｒの各要素に乗算してＪｒの行列を行と列とにＮ個展開すると、２進展開したＪｒは次の式（１５）のとおりとなる。
【００５３】
【数１５】
2020067811.tif 000017
【００５４】
ｈも同様に、重みδｗをＪｒの各要素に乗算してベクトルを行方向にＮ個展開し、これをｈ’とすると、ｈ’は次の式（１６）のとおりとなる。
【００５５】
【数１６】
2020067811.tif 000018
【００５６】
ピクセル（ｉ）については、Ｓの元となるＷは不変となる。Ｊ’、ｈ’については、各要素に対してｉ個分の要素を足し合わせたものとなり、次の式（１７）、（１８）のとおりである。なお、Ｓのビット数は、全部でＮ＊ｄ個となる。
【００５７】
【数１７】
2020067811.tif 000019
【００５８】
【数１８】
2020067811.tif 000020
【００５９】
最適化装置２００は、Ｓ３０４において、式（１７）、（１８）のＪ、ｈを結合係数（Ｑ）とし、Ｚ＝［αｗ，βｗ，δｗ］とともに出力する。なお、実際には、記号となっている箇所にはすべて数値が入る。
【００６０】
図１０に戻り、最適化装置２００は、結合係数（Ｑ）をもとにイジングマシンでシミュレーテッド・アニーリング（ＳＡ）による処理を行う（２０３）。シミュレーテッド・アニーリングでは、イジング形式に変換したコスト関数Ｆ’（Ｗ）の最小値が得られる各ビットの状態の組み合わせ（基底状態）を示すＷ’の数値を探索する。
【００６１】
次いで、最適化装置２００は、イジングマシンのシミュレーテッド・アニーリングにより得られたＷ’の数値に対して逆２進展開を行う（２０４）。これにより、最適化装置２００は、Ｗ’の数値に対応する、連続値であるＷを求める。具体的には、最適化装置２００は、２進展開の係数Ｚ＝［αｗ，βｗ，δｗ］を用いて、イジングスピン（Ｗ’の数値）から次の式（１９）によりＷｎ（連続値）を得る。
【００６２】
【数１９】
2020067811.tif 000021
【００６３】
ここで、Ｆ’（Ｗ）から得られた結果であるＷ’の値は正しい。しかしながら、不連続なｍａｘ関数となる項（Ｆ１（Ｗ））を２乗してｑ−ｌｏｓｓ関数に置き換えていることから、Ｆ（Ｗ）から見るとＷ’は近似値であり、正しい値とは言えない。
【００６４】
図１４は、ｑ−ｌｏｓｓ関数による変換例を説明する説明図である。具体的には、図１４に示すように、グラフＧ２０１は、不連続なｍａｘ関数となる項であるＦ１（ｍ）＝ｍａｘ（０，ｍ）の一例である。ｑ−ｌｏｓｓ関数を用いた置き換えでは、不連続なｍａｘ関数となる項Ｆ１（Ｗ）を二乗する。これにより、Ｆ１（ｍ）＝ｍａｘ（０，ｍ）２の一例であるグラフＧ２０２のように、関数の形（グラフの形状）が変わってしまう。
【００６５】
例えば、放射線照射計画の策定においては、照射するビームに対し、グラフＧ２０１に示す直線性に対応する正しい制約をつけることが重要である。しかしながら、グラフＧ２０２のように、変換後に直線性が保持されていない場合は、正しい制約をつけることが困難となり、ビームに対する制約が少なくなる、または過剰になることがある。
【００６６】
図１５は、不連続関数を例示する説明図である。図１５に示すように、イジングモデルに置き換えた際に関数の形が変わってしまう不連続関数としては、グラフＧ２０３の−ｍｉｎ（０，ｍ）や、グラフＧ２０４のｍａｘ（ｍ，−ｍ）などもある。例えば、グラフＧ２０３の−ｍｉｎ（０，ｍ）では、ｍを−ｍに変数変換すると、ｍａｘ（０，ｍ）となるため、ｍａｘ（０，ｍ）と同様に扱うことができる。
【００６７】
１つの側面では、より精度のよい最適解を求めることを可能とする最適化装置、最適化方法および最適化プログラムを提供することを目的とする。
【課題を解決するための手段】
【００６８】
１つの案では、最適化装置は、変換部と、２進展開部と、イジング処理部と、逆２進展開部とを有する。変換部は、入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換する。２進展開部は、第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする。イジング処理部は、第２のコスト関数の係数に関するイジング形式の値について、基底状態を示す値の探索を行う。逆２進展開部は、探索した基底状態を示す値について逆２進展開を行い、第１のコスト関数を最小とする値を求める。
【発明の効果】
【００６９】
本発明の１実施態様によれば、より精度のよい最適解を求めることができる。
【図面の簡単な説明】
【００７０】
【図１】図１は、実施形態にかかる最適化装置の機能構成例を示すブロック図である。
【図２】図２は、実施形態にかかる最適化装置の動作例を示すフローチャートである。
【図３】図３は、変換処理を例示するフローチャートである。
【図４】図４は、Ｌｅｇｅｎｄｒｅ変換の一例を説明する説明図である。
【図５】図５は、Ｗｏｌｆｅ双対定理を用いた変換の一例を説明する説明図である。
【図６】図６は、不連続関数について本実施形態と従来との比較例を示す説明図である。
【図７】図７は、実験結果の一例を説明する説明図である。
【図８】図８は、実施形態にかかる最適化装置のハードウェア構成例を示すブロック図である。
【図９】図９は、イジングマシンの一例を説明する説明図である。
【図１０】図１０は、従来の最適化装置の機能構成例を示すブロック図である。
【図１１】図１１は、放射線照射計画の策定例を説明する説明図である。
【図１２】図１２は、イジングモデルの式への変換の従来例を示すフローチャートである。
【図１３】図１３は、２進展開処理の一例を示すフローチャートである。
【図１４】図１４は、ｑ−ｌｏｓｓ関数による変換例を説明する説明図である。
【図１５】図１５は、不連続関数を例示する説明図である。
【発明を実施するための形態】
【００７１】
以下、図面を参照して、実施形態にかかる最適化装置、最適化方法および最適化プログラムを説明する。実施形態において同一の機能を有する構成には同一の符号を付し、重複する説明は省略する。なお、以下の実施形態で説明する最適化装置、最適化方法および最適化プログラムは、一例を示すに過ぎず、実施形態を限定するものではない。また、以下の各実施形態は、矛盾しない範囲内で適宜組みあわせてもよい。
【００７２】
図１は、実施形態にかかる最適化装置の機能構成例を示すブロック図である。図１に示すように、最適化装置１は、変換部１１、２進展開部１２、イジングマシン処理部１３および逆２進展開部１４を有する。
【００７３】
変換部１１は、最適化問題に関するコスト関数Ｆ（Ｗ）の入力を受け付け、Ｌｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いてイジングモデルの式（Ｆ’（Ｗ））への変換を行う処理部である。なお、本実施形態において、変換部１１に入力されるコスト関数Ｆ（Ｗ）は、一例として、図１１に例示した放射線照射計画に関するものとする。すなわち、入力されるコスト関数Ｆ（Ｗ）は、前述した式（３）のとおりである。
【００７４】
２進展開部１２は、変換したＦ’（Ｗ）が連続値のままであることから、イジングマシンで処理するために、Ｗを０，１に２進展開する処理部である。具体的には、従来の２進展開（２０２）と同様、２進展開部１２は、変換部１１が変換したコスト関数（Ｆ’（Ｗ））の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする。
【００７５】
イジングマシン処理部１３は、結合係数（Ｑ）をもとにイジングマシンでシミュレーテッド・アニーリング（ＳＡ）による処理を行い、基底状態を示す値の探索を行う処理部である。
【００７６】
逆２進展開部１４は、イジングマシン処理部１３が探索した基底状態を示す値、すなわちイジングマシンのシミュレーテッド・アニーリングにより得られたＷ’の数値に対して逆２進展開を行う。具体的には、逆２進展開部１４は、従来の逆２進展開（２０４）と同様の処理を行い、Ｗ’の数値に対応する連続値（Ｗ）を得る。これにより、逆２進展開部１４は、入力されたコスト関数Ｆ（Ｗ）を最小とする値を得る。
【００７７】
ここで、最適化装置１の各部の処理について詳細に説明する。図２は、実施形態にかかる最適化装置１の動作例を示すフローチャートである。
【００７８】
図２に示すように、処理が開始されると、変換部１１は、不連続項（Ｆ１（Ｗ））を持つコスト関数Ｆ（Ｗ）（式（３）参照）の入力を受け付ける（Ｓ１）。
【００７９】
ここで、不連続項（Ｆ１（Ｗ））は、式（３）に示すコスト関数Ｆ（Ｗ）における第２項（Σｉβｉｍａｘ（０，ｄｉ（ｗ）−Ｃｉ））である。この不連続項については、コスト関数Ｆ（Ｗ）の入力時にユーザーより指定されてもよいし、公知の数式処理システムによりコスト関数Ｆ（Ｗ）より検出してもよい。
【００８０】
次いで、変換部１１は、コスト関数Ｆ（Ｗ）の不連続項となる項（Ｆ１（Ｗ））について、Ｌｅｇｅｎｄｒｅ変換、ラグランジュ関数、およびＷｏｌｆｅ双対定理を用いてイジング形式に変換する（Ｓ２）。
【００８１】
ここで、Ｓ２の変換処理について詳細に説明する。図３は、変換処理を例示するフローチャートである。
【００８２】
図３に示すように、処理が開始されると、変換部１１は、コスト関数Ｆ（Ｗ）（式（３）参照）の入力を受け付ける（Ｓ１１）。次いで、変換部１１は、コスト関数Ｆ（Ｗ）において不連続項となる不連続関数Ｆ１（Ｗ）＝Σｉβｉｍａｘ（０，ｄｉ（ｗ）−Ｃｉ）について、ルジャンドル（Ｌｅｇｅｎｄｒｅ）変換を２回行う。これにより、変換部１１は、下に凸な関数Ｆ１の変数を変えて二次形式での式ｆ１（ｗ）を求める（Ｓ１２）。
【００８３】
不連続関数について、ルジャンドル変換した場合、傾きと切片の値で指定する接線の集合で同等に表現するので、変換しても性質は変わらない。例えば、関数ｆのルジャンドル変換ｆ＊（ｐ）は、次の式（２０）のとおりである。
【００８４】
【数２０】
2020067811.tif 000022
【００８５】
変換部１１は、不連続関数Ｆ１（Ｗ）のβｉｍａｘ（０，ｄｉ（ｗ）−Ｃｉ）に対してルジャンドル変換を２回行う。これは、ルジャンドル変換の正変換と逆変換は同じであるので、元に戻す操作に相当する。
【００８６】
例えば、ｆ１（ｗ）＝ｍａｘ（０，ｄｉ（ｗ）−Ｃｉ）＝−ｍｉｎ（０，Ｃｉ−ｄｉ（ｗ））とする。このｆ１にルジャンドル変換を２回かけると、ｍ＝ｄｉ（ｗ）−Ｃｉとして、ｆ１（ｍ）＝−ｍｉｎ（０，ｍ）と置くと次の式（２１）のとおりとなる。
【００８７】
【数２１】
2020067811.tif 000023
【００８８】
図４は、Ｌｅｇｅｎｄｒｅ変換の一例を説明する説明図である。図４において、グラフＧ１は、Ｌｅｇｅｎｄｒｅ変換前の不連続関数Ｆ１のグラフである。グラフＧ２は、不連続関数Ｆ１に対してＬｅｇｅｎｄｒｅ変換を２回かけた後のグラフである。なお、ｔは、Ｌｅｇｅｎｄｒｅ変換を２回かけた結果追加された変数である。グラフＧ２のｆ１（ｗ）の−（ｄｉ（ｗ）−Ｃｉ）を最小化（ｍｉｎ）するｔの条件は、ｄｉ（ｗ）＜Ｃｉでｔ＝０，ｄｉ（ｗ）＞Ｃｉでｔ＝−１となる。図４に示すように、不連続関数について、Ｌｅｇｅｎｄｒｅ変換を２回かけてもグラフの性質は変わらない。
【００８９】
ルジャンドル変換後の関数（ｆ１）のままでは、関数ｆ１のｍｉｎの前に−（マイナス）符号がついており、ｍとｔの最小化問題として取り扱うことができない。
【００９０】
したがって、変換部１１は、ルジャンドル変換後の関数（ｆ１）について、Ｗｏｌｆｅの双対定理により双対問題を求める（Ｓ１３）。すなわち、双対問題が解ければ、元の問題も解ける性質を利用する。
【００９１】
例えば、関数（ｆ１）について、双対問題を求めると次の式（２２）のとおりとなる。ここで、ｚ１、ｚ２はラグランジュ未定乗数である。ｍは、ｍ＝ｄｉ（ｗ）−Ｃｉである。
【００９２】
【数２２】
2020067811.tif 000024
【００９３】
また、式（２２）について、符号を逆にすると式（２３）のとおりとなる。
【００９４】
【数２３】
2020067811.tif 000025
【００９５】
図５は、Ｗｏｌｆｅ双対定理を用いた変換の一例を説明する説明図である。図５において、グラフＧ２は、図４と同様、不連続関数Ｆ１に対してＬｅｇｅｎｄｒｅ変換を２回かけた後のグラフである。グラフＧ３は、Ｗｏｌｆｅの双対定理により双対問題を求めて変換した式のグラフである。
【００９６】
図５に示すように、グラフＧ３の式では、Ｌｅｇｅｎｄｒｅ変換を２回かけた際に追加されたｔに、ｚ１，ｚ２が追加される。グラフＧ３において、ｆ１（ｗ）を最小化（ｍｉｎ）するｚ１，ｚ２の条件は、グラフの形に添わせるためｄｉ（ｗ）＜Ｃｉではｚ１＝０となる。また、ｄｉ（ｗ）＞Ｃｉでは、ｚ１，ｚ２は不定、すなわちグラフの形に関与しないということになる。
【００９７】
図３に戻り、Ｓ１３に次いで、変換部１１は、Ｗｏｌｆｅの双対の制約条件（双対問題を求めて変換した式）をペナルティ項（Ｍ（−ｄｉ（ｗ）−Ｃｉ）−ｚ１＋ｚ２）２として２体イジング形式で表現する（Ｓ１４）。具体的には、次の式（２４）のとおりとなる。なお、Ｍは、十分に大きな数値（例えば１００など）とする。
【００９８】
【数２４】
2020067811.tif 000026
【００９９】
次いで、変換部１１は、コスト関数（Ｆ（Ｗ））について、不連続関数（Ｆ１（Ｗ））に代わり、２体イジング形式で表現したペナルティ項を盛り込んだ関数（Ｆ’（Ｗ））全体と値の範囲を出力する。
【０１００】
具体的には、出力するＦ’（ｗ）は、次の式（２５）のとおりである。また、出力するｚ１，ｚ２，ｔの値の範囲は、ｚ１≧０，ｚ２≧０，−１≦ｔ≦０となる。
【０１０１】
【数２５】
2020067811.tif 000027
【０１０２】
図２に戻り、Ｓ２において、イジング形式に変換したイジングスピンは、実数値である。したがって、Ｓ２に次いで、２進展開部１２は、イジング形式の実数値のスピンに対して２進展開を行い、２値のスピンでのイジング形式を求める（Ｓ３）。具体的には、２進展開部１２は、従来の２進展開（２０２）と同様、２進展開の精度（ｄｗ）、値の範囲（α、β）の入力を受け付け、Ｗの２進展開を行う。次いで、２進展開部１２は、２進展開を行ったＷ’に関する結合整数（ｑ）およびＷ’を元の実数に戻すための係数（Ｚ）を出力する。
【０１０３】
次いで、イジングマシン処理部１３は、イジング形式を用いてイジングマシンによる処理を行う（Ｓ４）。具体的には、イジングマシン処理部１３は、結合係数（Ｑ）をもとにイジングマシンでシミュレーテッド・アニーリング（ＳＡ）による処理を行い、基底状態を示す各スピンの値を求める。
【０１０４】
次いで、逆２進展開部１４は、イジングマシンによる処理結果のスピンの値（基底状態を示す各スピンの値）を用いて、２進展開前の係数の値を求める（Ｓ５）。具体的には、逆２進展開部１４は、従来の逆２進展開（２０４）と同様の処理を行い、イジングマシンより得られたＷ’の数値に対して逆２進展開を行う。これにより、最適化装置１は、入力されたコスト関数Ｆ（Ｗ）を最小とする値（実数値）を得る。
【０１０５】
以上のように、最適化装置１は、変換部１１と、２進展開部１２と、イジングマシン処理部１３と、逆２進展開部１４とを有する。変換部１１は、入力された第１のコスト関数（Ｆ（Ｗ））をラグランジュ関数とＷｏｌｆｅの双対定理とを用いて第２のコスト関数（Ｆ’（Ｗ））に変換する。２進展開部１２は、第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする結合係数（Ｑ）を求める。イジングマシン処理部１３は、第２のコスト関数（Ｆ’（Ｗ））に関する結合係数について、基底状態を示す値の探索を行う。逆２進展開部１４は、探索した基底状態を示す値について逆２進展開を行い、第１のコスト関数（Ｆ（Ｗ））を最小とする値を求める。
【０１０６】
これにより、最適化装置１は、第１のコスト関数（Ｆ（Ｗ））に含まれる絶対値制約をイジングモデルで記述する際に、関数の形状（グラフ形状）を変えることなく、絶対値制約のまま扱うことができる。したがって、最適化装置１では、より精度のよい最適解を求めることができる。
【０１０７】
図６は、不連続関数について本実施形態と従来との比較例を示す説明図である。図６において、グラフＧ４は、不連続関数Ｆ１（ｍ）を変換部１１によりイジングモデルで記述した場合のグラフである。グラフＧ４は、不連続関数Ｆ１（ｍ）をｑ−ｌｏｓｓ関数を用いた従来の方法によりイジングモデルで記述した場合のグラフである。
【０１０８】
図６に示すように、グラフＧ４では、ｍａｘ関数に対応するグラフ形状のまま（絶対値制約のまま）扱うことができる。これに対し、従来の方法では、変換後のグラフ形状において直線性が保持されておらず、制約が少なくなる領域Ｇ５１や、制約が過剰になる領域Ｇ５２が生じる。この結果、例えば、放射線照射計画の策定においては、ビームに対する制約が少なくなる、または過剰になる場合がある。
【０１０９】
図７は、実験結果の一例を説明する説明図である。この図７に示す実験結果は、腫瘍（病Ｐ（１）〜（４））と、正常組織（正Ｐ（１））とをピクセル配置Ｐ１として配置した場合の放射線照射計画を求めたものである。なお、ビームは時計回りに２０箇所あり、照射量（Ｍ）は、行方向にＩ、列方向にＮとして図示のとおりである。また、Ｐ（１〜５）の値、およびＣｉの値は図示のとおりとする。また、得られた強度ＷとＭとを乗算した結果がＣｉを超えなければ、ｑ−ｌｏｓｓ関数を用いた二乗制約でもかまわないものとする（ここではＣｉ＝Ｐ＋０．０１）。
【０１１０】
図７において、ＭＷａｂｓは、図示の条件において、本実施形態の最適化装置１より（絶対値制約で）得られたＷとＭとを乗算した結果である。また、ＭＷｓｑｕａｒｅは、従来の最適化装置２００より（二乗制約で）得られたＷとＭとを乗算した結果である。図７に示すように、ＭＷｓｑｕａｒｅでは、１つ目と３つ目の値がＣｉより大きく、過剰になっている。これに対し、ＭＷａｂｓでは、Ｃｉを超える値がない。したがって、本実施形態の最適化装置１による放射線照射計画では、ビームに対する制約が有効に働いていることがわかる。
【０１１１】
また、変換部１１は、第１のコスト関数（Ｆ（Ｗ））に含まれる不連続関数（Ｆ１（Ｗ））について、ルジャンドル変換を行う。そして、変換部１１は、ルジャンドル変換後の式をもとにＷｏｌｆｅの双対定理による双対問題を求めて得られた関数を、不連続関数に代えて第１のコスト関数に盛り込むことで第２のコスト関数（Ｆ’（Ｗ））に変換する。これにより、最適化装置１は、第１のコスト関数（Ｆ（Ｗ））に含まれる不連続関数の形状（グラフ形状）を変えることなく、第１のコスト関数（Ｆ（Ｗ））から第２のコスト関数（Ｆ’（Ｗ））に変換することができる。
【０１１２】
なお、図示した各装置の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。
【０１１３】
最適化装置１で行われる各種処理機能は、ＣＰＵ（Central  Processing  Unit）（またはＭＰＵ（Micro-Processing  Unit）、ＭＣＵ（Micro  Controller  Unit）等のマイクロ・コンピュータ）上で、その全部または任意の一部を実行するようにしてもよい。また、各種処理機能は、ＣＰＵ（またはＭＰＵ、ＭＣＵ等のマイクロ・コンピュータ）で解析実行されるプログラム上、またはワイヤードロジックによるハードウェア上で、その全部または任意の一部を実行するようにしてもよいことは言うまでもない。また、最適化装置１で行われる各種処理機能は、クラウドコンピューティングにより、複数のコンピュータが協働して実行してもよい。
【０１１４】
ところで、上記の実施形態で説明した各種の処理は、予め用意されたプログラムをコンピュータで実行することで実現できる。そこで、以下では、上記の実施例と同様の機能を有するプログラムを実行するコンピュータ（ハードウエア）の一例を説明する。図８は、実施形態にかかる最適化装置１のハードウエア構成例を説明する説明図である。
【０１１５】
図８に示すように、最適化装置１は、入出力装置１００と、記憶装置１０１と、ＲＡＭ１０２（ＲＡＭ：Random  Access  Memory）と、ＲＯＭ１０３（ＲＯＭ：Read  Only  Memory）と、ＣＰＵ１０４と、イジングマシン１０５とを有する。最適化装置１内の各部（１００〜１０５）は、バス１０６に接続される。
【０１１６】
入出力装置１００は、ユーザーからの入力操作を受け付けるキーボード、マウスなどの入力装置、および、各種処理結果を出力するディスプレイなどの出力装置である。記憶装置１０１は、例えばハードディスク装置などであり、上記の実施形態で説明した各種の処理を実行するためのプログラム１０１ａが記憶される。また、記憶装置１０１には、プログラム１０１ａが参照する各種データが記憶される。ＲＡＭ１０２は、ＣＰＵ１０４がプログラム１０１ａを読み出して各種の処理を実行する際に用いられ、各種情報を一時記憶する。ＲＯＭ１０３は、例えば最適化装置１の起動時に実行されるブートプログラム等を記憶する不揮発性メモリである。
【０１１７】
ＣＰＵ１０４は、記憶装置１０１に記憶されたプログラム１０１ａを読み出してＲＡＭ１０２に展開して実行することで、変換部１１、２進展開部１２、イジングマシン処理部１３および逆２進展開部１４に関する各種の処理を行う。なお、プログラム１０１ａは、記憶装置１０１に記憶されていなくてもよい。例えば、最適化装置１が読み取り可能な記憶媒体に記憶されたプログラム１０１ａを、最適化装置１が読み出して実行するようにしてもよい。最適化装置１が読み取り可能な記憶媒体は、例えば、ＣＤ−ＲＯＭやＤＶＤディスク、ＵＳＢ（Universal  Serial  Bus）メモリ等の可搬型記録媒体、フラッシュメモリ等の半導体メモリ、ハードディスクドライブ等が対応する。また、公衆回線、インターネット、ＬＡＮ等に接続された装置にこのプログラム１０１ａを記憶させておき、最適化装置１がこれらからプログラム１０１ａを読み出して実行するようにしてもよい。
【０１１８】
イジングマシン１０５は、イジングマシン処理部１３の制御のもと、シミュレーテッド・アニーリング（ＳＡ）による処理を行い、イジング形式に変換したコスト関数の最小値が得られる各ビットの状態の組み合わせ（基底状態）を探索する装置である。
【０１１９】
図９は、イジングマシン１０５の一例を説明する説明図である。図９に示すように、イジングマシン１０５は、コントローラ１１０と、ハードウェア最適化エンジン１１１と、複数のニューロン１１２とを有する。
【０１２０】
コントローラ１１０は、結合係数に関する重みＷｉｊをハードウェア最適化エンジン１１１に渡してシミュレーテッド・アニーリング（ＳＡ）処理を行わせ、各ビットの状態を得る。ハードウェア最適化エンジン１１１は、コントローラ１１０からの重みＷｉｊをニューロン１１２各々に設定し、ニューロン１１２各々よりＳＡ処理後の出力を得る。次いで、ハードウェア最適化エンジン１１１は、ニューロン１１２各々より得られた出力をコントローラ１１０に出力する。
【０１２１】
複数のニューロン１１２は、ニューラルネットワークを用いてイジングモデルをモデル化したものである。具体的には、複数のニューロン１１２のそれぞれが、他のビットの状態と、他のビットと自身のビットとの結合の強さを示す重みＷｉｊ（結合係数）とに応じて０または１を出力するニューロンとして機能する。
【０１２２】
具体的には、ニューロン１１２各々は、乗算器１１３と、加算器１１４と、コンパレータ１１６とを有する。乗算器１１３のそれぞれは、ハードウェア最適化エンジン１１１から出力されるｎ個のニューロン１１２の状態（ｘ１，ｘ２，…，ｘｎ）と、重みＷｉ１，…，Ｗｉ，ｎとを掛け合わせる。加算器１１４は、乗算器１１３から出力される値を積算し、ノイズ１１５を加える。コンパレータ１１６は、加算器１１４が積算してノイズ１１５を加えた値と、閾値との比較結果を出力する。ニューロン１１２は、コンパレータ１１６より得られた出力値をもとにニューロン１１２のパラメータを順次更新して処理を繰り返すことで、基底状態を探索する。
【０１２３】
以上の実施形態に関し、さらに以下の付記を開示する。
【０１２４】
（付記１）入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換する変換部と、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とする２進展開部と、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行うイジング処理部と、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める逆２進展開部と、
を有することを特徴とする最適化装置。
【０１２５】
（付記２）前記変換部は、前記第１のコスト関数に含まれる不連続関数について、ルジャンドル変換を行い、当該ルジャンドル変換後の式をもとにＷｏｌｆｅの双対定理による双対問題を求めて得られた関数を、前記不連続関数に代えて前記第１のコスト関数に盛り込むことで、前記第２のコスト関数に変換する、
ことを特徴とする付記１に記載の最適化装置。
【０１２６】
（付記３）入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換し、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とし、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行い、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める、
処理をコンピュータが実行することを特徴とする最適化方法。
【０１２７】
（付記４）前記変換する処理は、前記第１のコスト関数に含まれる不連続関数について、ルジャンドル変換を行い、当該ルジャンドル変換後の式をもとにＷｏｌｆｅの双対定理による双対問題を求めて得られた関数を、前記不連続関数に代えて前記第１のコスト関数に盛り込むことで、前記第２のコスト関数に変換する、
ことを特徴とする付記３に記載の最適化方法。
【０１２８】
（付記５）入力された第１のコスト関数をＬｅｇｅｎｄｒｅ変換、ラグランジュ関数、および、Ｗｏｌｆｅ双対定理を用いて第２のコスト関数に変換し、
前記第２のコスト関数の係数に関する値を、２進展開を行って２値のスピンでのイジング形式とし、
前記第２のコスト関数の係数に関する前記イジング形式の値について、基底状態を示す値の探索を行い、
探索した前記基底状態を示す値について逆２進展開を行い、前記第１のコスト関数を最小とする値を求める、
処理をコンピュータに実行させることを特徴とする最適化プログラム。
【０１２９】
（付記６）前記変換する処理は、前記第１のコスト関数に含まれる不連続関数について、ルジャンドル変換を行い、当該ルジャンドル変換後の式をもとにＷｏｌｆｅの双対定理による双対問題を求めて得られた関数を、前記不連続関数に代えて前記第１のコスト関数に盛り込むことで、前記第２のコスト関数に変換する、
ことを特徴とする付記５に記載の最適化プログラム。
【符号の説明】
【０１３０】
１…最適化装置
１１…変換部
１２…２進展開部
１３…イジングマシン処理部
１４…逆２進展開部
１００…入出力装置
１０１…記憶装置
１０１ａ…プログラム
１０２…ＲＡＭ
１０３…ＲＯＭ
１０４…ＣＰＵ
１０５…イジングマシン
１１０…コントローラ
１１１…ハードウェア最適化エンジン
１１２…ニューロン
１１３…乗算器
１１４…加算器
１１５…ノイズ
１１６…コンパレータ
２００…最適化装置
３０１…腫瘍
３０２…正常組織
３０３…マルチリーフコリメータ
Ｍ…影響度
Ｐ１…ピクセル配置
Ｇ１〜Ｇ５、Ｇ２０１〜Ｇ２０４…グラフ
Ｇ５１、Ｇ５２…領域
【図１】
2020067811.tif 000028
【図２】
2020067811.tif 000029
【図３】
2020067811.tif 000030
【図４】
2020067811.tif 000031
【図５】
2020067811.tif 000032
【図６】
2020067811.tif 000033
【図７】
2020067811.tif 000034
【図８】
2020067811.tif 000035
【図９】
2020067811.tif 000036
【図１０】
2020067811.tif 000037
【図１１】
2020067811.tif 000038
【図１２】
2020067811.tif 000039
【図１３】
2020067811.tif 000040
【図１４】
2020067811.tif 000041
【図１５】
2020067811.tif 000042