JP2003030172A

JP2003030172A - 拘束条件付き最適化方法及びプログラム

Info

Publication number: JP2003030172A
Application number: JP2001210796A
Authority: JP
Inventors: Akira Asai; 朗浅井; Shigeki Matsutani; 茂樹松谷
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2001-07-11
Filing date: 2001-07-11
Publication date: 2003-01-31
Also published as: US7197486B2; US20030065689A1

Abstract

(57)【要約】【課題】拘束条件付き最適化問題の解を簡単な計算手
順で求める。【解決手段】等式で与えられる拘束条件の下での最適
化関数の最小値を求める拘束条件付き最適化方法におい
て、条件式で定義される集合を有限次元実ベクトル空間
内のリーマン多様体と見なし、実ベクトル空間内の初期
位置x₀から、拘束条件を満たす当該リーマン多様体へ接
近し(S2)、前記リーマン多様体上の接ベクトルに対する
測地線方程式に関する指数関数写像を有限次数で打ち切
り、該打ち切りによって得られた近似測地線を実ベクト
ル空間内の１次元軌道として生成し、前記リーマン多様
体上の接ベクトルと生成される前記軌道に対して当該接
ベクトルの平行移動に関する指数関数写像の有限次数近
似によって近似平行移動を行う(S8)。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、等式で与えられる
拘束条件の下での最適化関数の最小値を求める最適化方
法、その装置、及びそのプログラムに関するものであ
る。

【０００２】

【従来の技術】等式で指定された拘束条件 S(r)＝0， i=1，・・・，m (1) の下で、実数値関数E(r)を極小化する問題を考える。こ
こでｒはn次元実ベクトル即ち、

【外２】

【０００３】であるとする。なお、実関数を極大化する
問題は、関数の符号を逆にすることにより、上述の問題
に帰着する。

【０００４】拘束条件のない場合、このような問題を解
くためのアルゴリズムとしては、適当な初期ベクトルｒ
⁽⁰⁾から出発し、勾配ベクトルF = ―∇E(r)に関する情
報を用いて、関数E(r)が極小になる点ｒを探索するもの
が一般的である。

【０００５】近年ではFletcher-Reeves法、Polak−Ribi
ere法等を含む共役勾配法も用いられる。これらのアル
ゴリズムの詳細は(「Computational Methods in Optimi
zation」 E.Polak著 1971年Academic Press発行)に書か
れている。これらの方法では前回探索した方向と垂直に
探索方向を選ぶために、探索の効率化を行え、関数E(r)
が２次近似できるところでは、Hestnes−Stiefelの共役
勾配法の代数的なよい性質を利用できる利点があり、近
年多く用いられている。

【０００６】拘束条件の加わる場合は、このような単純
なアルゴリズムでは巧くいかないので、種々の工夫がな
されてきた。代表的なものとしては、ペナルティ法（罰
金法）と乗数法（ラグランジュの未定乗数法）が良く知
られている。これらの方法等に関しては、例えば「岩波
講座応用数学、最適化法、藤田宏、今野浩、田邉國士
著、1994年岩波書店発行」に詳しく書かれている。

【０００７】ペナルティ法は、拘束条件を表現するため
に、拘束条件が成り立つときには０となり、成り立たな
ければ非常に大きな数になるようなペナルティ関数（罰
金関数）P(r)を導入し、新たな関数 Ω(r)＝E(r)＋P(r) (3) を極小化するものである。ペナルティ関数の具体的な形
としては、

【外３】

【０００８】などが考えられる。ここにω_kは正の数
で、極小値探索の過程で、適当な値に調整される。

【０００９】乗数法は、未定乗数λ_kを導入し、ｒとλ_k
の新たな関数

【外４】

【００１０】を停留化するものである。拘束条件は、λ
_kに対するΩ(r，λ_i)の停留条件

【外５】

【００１１】として実現される。

【００１２】また、これらの方法とは別に、考えている
空間の中にベクトル場を考え、最適点をその漸近解とす
るような力学系を用いたアルゴリズムも研究されている
（上述の岩波講座「最適化法」に述べられている）。

【００１３】更に、最近、Smithらは、許容集合すなわ
ち拘束条件をみたす点の集合をリ−マン多様体と捉え、
その上でのニュートン法および共役勾配法を考えるとい
うアルゴリズムを提案している。（S.T. Smith， AMS，
Fileds Institute Communications vol.3， 1994， p
p.113-136，「Optimaization Techniques on Riemanni
an Manifolds」）また、許容集合が Grassmann 多様体
またはStiefel多様体となる場合についてはより詳しい
アルゴリズムを提案している。(A. Edelman， T.A. Ari
as， S.T.Smith， SIAM J. Matrix Anal. Appl. vol. 2
0， 1998， pp.303-353、「The Geometry of Algorithm
s with Orthogonality Constraints」)アルゴリズムは
拘束条件のない場合と基本的に同様になるので、他の方
法に現れるような追加変数や調整パラメータが不要であ
るという特長がある。この手法は、拘束条件付最適化問
題に対する新しい視点を提供するもので、固有値問題な
どへの幅広い適用が期待できる。

【００１４】

【発明が解決しようとする課題】従来の方法で述べたよ
うに、等式拘束条件のもとでの最適化問題のアルゴリズ
ムとしては、乗数法と罰金法が広くもちいられている。
前者は鞍点探索問題となり、極値探索に比べて計算アル
ゴリズムが面倒になる。また、後者は罰金のパラメータ
選択を問題に応じて工夫しなければならないという問題
点がある。

【００１５】また、上述した力学系を用いたアルゴリズ
ムでは、適当な力学系を設定するための工夫を必要とす
る。

【００１６】また、上記Smithらのアルゴリズムにおい
ては、最適点への探索過程において、点を許容集合上の
測地線に沿って移動させる。Smithらは、Grassmann 多
様体とStiefel多様体の場合について、測地線方程式の
解を具体的に構成している。原理的には、この解に従っ
て点を移動させると、点は許容集合内に留まるはずだ
が、現実的には数値計算上の誤差のため、点が許容集合
から離れるという問題がある。また、許容集合がこれら
の多様体でない場合については、測地線方程式の具体的
な解の表式が与えられていない。

【００１７】

【課題を解決するための手段】上記課題を解決するため
に、本発明によれば、等式で与えられる拘束条件の下で
の最適化関数の最小値を求める拘束条件付き最適化方法
に、条件式で定義される集合を有限次元実ベクトル空間
内のリーマン多様体と見なし、実ベクトル空間内の初期
位置から、拘束条件を満たす当該リーマン多様体へ接近
する接近工程と、前記リーマン多様体上の接ベクトルに
対する測地線方程式に関する指数関数写像を有限次数で
打ち切り、該打ち切りによって得られた近似測地線を実
ベクトル空間内の１次元軌道として生成する軌道生成工
程と、前記リーマン多様体上の接ベクトルと前記軌道生
成工程によって生成される軌道に対して当該接ベクトル
の平行移動に関する指数関数写像の有限次数近似によっ
て近似平行移動を行う近似平行移動工程とを備える。

【００１８】

【発明の実施の形態】ここで、上述した等式で与えられ
る拘束条件の下での最適化関数の最小値を求めるという
問題を数学的に定式化する。

【００１９】n次元実ベクトル空間Ｒⁿを考え、(e_i|i=
1，・・・，n)をその基底とする。(e_i|i=1，・・・，n)は必ず
しも直交していないとし、内積< ， >に対して b_ij=<e_i， e_j> i，j=1，・・・，n (7) の値を取るとしている。実際、有限要素法で対象とする
最適化問題においては、基底は直行化されていない。更
に、m個の拘束条件(1≦m≦n)を表わすＣ^∞級関数、 S^k：Ｒⁿ→Ｒ (k=1，・・・，m) (8) および下に有界なＣ^∞級関数（目的関数）Ｅ：Ｒⁿ→Ｒ (9) が与えれらているとする。拘束条件は上記の関数を使っ
てS^k(x)=0， k=1，・・・，mとして与えられているとす
る。、拘束条件を満たすものを許容集合と称しＭ={x∈Ω|S^k(x)=0， for k=1，・・・，m } (10) と記す。このとき、ＭのＲⁿ内の近傍において正則条件 rankJ(x)=m (11) を満たすものとする。ここにJ(x)はヤコビ行列で、その
ki成分は、

【外６】

【００２０】で定義される。

【００２１】ユークリッド空間の部分空間Ｒⁿはbを計量
とする平坦なリ−マン多様体とみなすことができる。拘
束条件の正則性から、ＭもＲⁿの計量をＭ上に制限する
ことによって得られるb|Ｍを計量としてＲⁿのm次元部分
リーマン多様体となる。従って、考えている拘束条件付
き最適化問題は、リ−マン多様体Ｍ上の最適化問題：
「Ｅ(x)を最小化するx∈Ｍを求めよ」として定式化され
る。

【００２２】次にＲⁿにおけるベクトル解析を利用し
て、多様体Ｍを記述することを示す。n-m次元リ−マン
多様体上の点を表現するための座標としては局所座標
(ξⁱ|i＝1，・・・，n-m)を用いる事ができるが、ここで
は、ＭがＲⁿの部分多様体であることを利用して、Ｒⁿの
位置ベクトルを用いてＭを表現する。この座標(xⁱ|i＝
1，・・・，n-m)をＭにおける局所座標と区別するため、こ
こでは外的座標（extrinsic coordinate）と呼ぶことに
する。数値計算の立場からは、外的座標を用いると、線
形代数に関するアルゴリズムの蓄積を用いて効率的なア
ルゴリズムを構築できるという利点がある。リ−マン多
様体としてのＭ及びその近傍の性質を外的座標を用いて
記述する。

【００２３】そこで、ここでは、そのための道具として
ＲⁿにおけるGradient、方向微分、Hessianを定義してお
く。これらの、Ｒⁿにおけるベクトル解析の表現を用い
ると外的座標の選択によらない形でＭ上の量を記述する
ことが可能となり、見通しがよくなる。

【００２４】なお、以下の記述において、外的座標に関
する添字に対しては、アインシュタインの規約を用い
る。つまり、同じ添え字が２回出て来た際には、その添
え字について1からnまでの和を取るものとする。また、
行列(b_ij|i，j＝1，・・・，n)の逆行列のij成分をb^ijと書
くことにする。

【００２５】以下にGradient、方向微分、Hessianを定
義する。

【００２６】xをＲⁿ内の点、fをＲⁿ上のＣ^∞級関数f：
Ｒⁿ →Ｒとするとき、任意のベクトルZ＝Zⁱe_iに対し
て、

【外７】

【００２７】を満たすベクトルYが一意的に存在する。
このベクトルを、fのxにおけるgradient（勾配）と呼
び、Grad_xfと書く。成分を用いて表示すると、

【外８】

【００２８】となる。計算の便宜上のために以下の量を
導入する。

【００２９】Ｆ^k(x):＝ Grad_xS^k， (k＝1，・・・，m) (15) Ｑ^kl(x):＝ <Ｆ^k(x)，Ｆ^l(x)>， (k，l＝1，・・・，m) Ｑ_kl(x):＝行列(Ｑ^kl(x)| k，l＝1，・・・，m)の逆行列のk，l成分但し、行列Ｑ_kl(x)の存在は、Ｑ^kl(x)の行列としての拘
束条件の正則性に依存するが、拘束条件の正則性は許容
集合Ｍ上では保証されておりＲⁿ内Ｍの近傍においても
保証されている。多くの物理的な例においては、関数S^k
(x)は充分穏やかな関数であり、Ｍ自身も充分穏やかな
形状を持つため、以下で考察する領域内においてはＱ_kl
(x)は存在しているとして良い。次に方向微分を導入する。Ｒⁿ内の点xにおいて、関数f
∈Ｃ^∞(Ω)の、ベクトルY=Yⁱe_iにそっての方向微分(Y・G
rad)_xfを、

【外９】

【００３０】と定義する。また、ベクトル場Z(x)＝Z
ⁱ(x)e_iの、ベクトルY＝Yⁱe_iにそっての方向微分(Y・Gra
d)_xZを、

【外１０】

【００３１】によって定義する。この方向微分の定義
は、外的座標の選び方によらないことは良く知られた事
実である。

【００３２】次にHessianを定義する。xをＲⁿ内の点、f
をＲⁿ上のＣ^∞級関数f:Ω→Ｒとするとき、ベクトルY=Y
ⁱe_i， Z=Zⁱe_iを引数とする実数値関数Hess_xfを

【外１１】

【００３３】によって定義する。この関数Hess_xfをfのx
におけるHessianと呼ぶ。

【００３４】次に上述のGradient、方向微分、Hessian
の定義を利用して、許容集合Ｍへの接近アルゴリズムを
述べる。xをＲⁿ内の点でＭの近傍の点であるとする。近
傍系においては、関数の族{S^k}をＭをゼロ点とする法線
方向のベクトル空間の成分と近似できる。つまり、Ｍ上
の点yをＭ上の点で、xに最も近い点とすると良い精度で

【外１２】

【００３５】と近似できる。但し、μ_kはそれぞれn次元
のベクトル量で成分表示するとμ_k =(μ_k ¹，μ_k ²，・・
・，μ_k ⁿ)となる。yをＭ上の点としたので、S^k(y)=0であ
る。仮定より

【外１３】

【００３６】であるので、

【外１４】

【００３７】テイラー展開の一次で近似すると、

【外１５】

【００３８】となる。式(15)の量を援用すると

【外１６】

【００３９】となる。近傍系の一次独立性を考慮すると <F^k(x)，μ_l>=δ_kl (21) となる。この式をμ_kについて解くと、式(15)より

【外１７】

【００４０】となる。従って、Ｒⁿ内の点xでＭの近傍の
点において、xのＭへの射影を行った際にその変化分の
ベクトルを

【外１８】

【００４１】と定義する。

【００４２】次に埋め込まれた多様体の性質として良く
知られた技法を用いて、許容集合Ｍのリ−マン幾何学的
性質を述べる。積分可能性の議論においてよく知られて
いるようにＭ上のベクトルYがＭの接ベクトルとなるた
めの必要十分条件は、 (dS^k)_x(Y)＝0 または<Grad_xS^k，Y> (24) である。Ｍ上の点ｘにおいてＲ^ｎにおけるベクトルYが
与えられたとき、ｘにおける任意のＭの接ベクトルＵ∈
Ｔ_xＭに対して <Ｕ，Ｖ>＝<Ｕ，Ｙ> (25) をみたすＶ∈Ｔ_xＭが一意的に存在し、

【外１９】

【００４３】と書くことができる。ベクトルYの接空間
Ｔ_xＭへの射影ＶをProj_xYと書く。

【００４４】Ｍにおける接ベクトル場Ｕ，Ｖ∈χ(Ｍ)が
与えられたとき、Ｕ方向へのＶの共変微分∇_UＶ∈χ
(Ｍ) は、

【外２０】

【００４５】によって与えられる。これを式(26)のProj
_xの定義を用いて書き直すと、

【外２１】

【００４６】となる。

【００４７】ここでＶがＭにおける接ベクトル場である
ことより得られる関係式<F^k，Ｖ>＝0にＵ・Gradを施す
ことにより、という関係式を得る。これを援用すると ∇_UＶ＝(Ｕ・Grad)_xＶ―Γ(Ｕ，Ｖ) (30) 但し、

【外２２】

【００４８】である。

【００４９】次に平行移動について述べる。Ｍ上のベク
トル場Ｖが滑らかな曲線

【外２３】

【００５０】に沿って平行であるための条件は、曲線上
の各点γ(t)において曲線の接線方向へのＶの共変微分
が0、すなわち

【外２４】

【００５１】である。ここで、

【外２５】

【００５２】であることを用いると、良く知られたＭ上
の曲線

【外２６】

【００５３】上のベクトル場V(t)の平行移動方程式、

【外２７】

【００５４】となる。

【００５５】初期条件V(0)=Ｇに対する平行移動方程式
の解は、次の式で与えられる： V(t)=Ｇ-tΓ(H，G)+o(t) (35) ここに、

【外２８】

【００５６】とした。Ｍ上の曲線

【外２９】

【００５７】が測地線であるとは、

【外３０】

【００５８】が平行移動方程式を満たすことである。従
って、測地線方程式は、

【外３１】

【００５９】となる。初期条件γ(0)=x₀に対する測地線
方程式の解は、次の式で与えられる：

【外３２】

【００６０】ここに、

【外３３】

【００６１】である。

【００６２】次にSmithのアルゴリズムについて概略を
述べる。そもそも一般的にリーマン多様体上において方
向という概念は巧く定義されたものではない。リーマン
多様体上のある点での接ベクトルを定義し、その点を通
るある閉曲線上で平行移動させた際、初期ベクトルと閉
曲線を回ってきたベクトルの方向は一般に異なる。その
変化分は採用した閉曲線内のスカラー曲率に依存する事
が良く知られている。

【００６３】その為、多様体上で共役勾配法を構成する
という試み自体、2次形式に対する共役勾配法のもつ良
い性質の全てが満たされることは期待できない。しかし
ながら、１）探索という操作をＭへの曲線（1次元部分多様体）
の埋め込みと考えると接線ベクトル

【外３４】

【００６４】に対する垂直方向というものは巧く定義さ
れる。よって、新たな探索方向を

【外３５】

【００６５】に垂直になるように補正することは、無駄
の無い計算方法と考えられ、意味を持つ。２）また、拘束条件の下で最小値に充分近づいた際でか
つ、許容多様体Ｍ上での最小値近傍ではＭに沿ってのエ
ネルギー関数が充分２次形式であると見なすことができ
るという場合に限られるが、その際にはHestenes，Lanc
zos，Stiefelの共役勾配法が実効的に働き、収束を早め
ることが期待できる。

【００６６】よって、エネルギー積分が２次形式の場
合、Hestenes， Lanczos， Stiefelと同様のアルゴリズ
ムになるように、多様体上の共役勾配法なるものを定義
することは自然な考えである。実際、物理現象をモデル
化した系の場合はエネルギー関数がそのような性質をも
つ場合が多いことが経験的に知られている。そこで、Sm
ithのアルゴリズムについて図７に従って述べる。

【００６７】Smithのアルゴリズム：ステップS71で初期
位置設定をx=x₀とする。

【００６８】ステップS72で、許容多様体Ｍ上に点が無
ければ、許容多様体Ｍ上の点へ移動させる。

【００６９】ステップS73で、エネルギーの収束を判断
する。一度目はこれがNoとなるようにしておく。

【００７０】ステップS74で、Ｇ_new:=Proj_xGrad_xE (38) を計算する。

【００７１】ステップS75で反復回数を判定する。初回
であれば、ステップS76に移行し、Ｇ:=Ｇ_new、Ｈ:=-Ｇ
_newとする。

【００７２】ステップS78で１）H方向の測地線上で、E(t)が最小となるtを求める。２）ベクトルGを測地線に沿ってtだけ平行移動させたも
のをτ₀(Ｇ)とする。３）ベクトルHを測地線に沿ってtだけ平行移動させたも
のをτ₀(Ｈ)とする。４）点xを測地線上をtだけ移動させる。

【００７３】そして、ステップS73に戻って、収束判定
を行う。エネルギーの変化量が所望の値より小さい場合
は収束したと判断して、終了する。収束していない場合
は、ステップS74でＧ_newを計算した後、ステップS75の
判定を行った後にステップS77の手順を踏む。

【００７４】

【外３６】

【００７５】

【外３７】

【００７６】

【外３８】

【００７７】その後再び、ステップS78へ移行する。上
記のステップを流れ図に従い、計算する。

【００７８】しかしながら、Smithのアルゴリズムで数
値計算を行うと、数値誤差のため、平行移動の結果得ら
れるτ₀(Ｇ)、τ₀(Ｈ)は接ベクトルではなくなり、ま
た、xは、Ｍ上の点ではなくなってしまう。そのため、
ステップS71で行った許容多様体Ｍへの接近をステップS
78においても毎回行わなければならない。そもそも、上
述したように、このアルゴリズムは一般の拘束条件やエ
ネルギー関数に対しては、共役勾配法として数学的に正
当化されるものではなく、平行移動等を正確に計算する
事に重大な意味があるものではない。

【００７９】他方、τ₀(Ｇ)、τ₀(Ｈ)等の操作は計算機
コストのかかる計算であり、計算時間の浪費を招いてい
る。そこで、本実施形態においては、この点を以下の２
次近似で置き換えることで代用する事を提案する。許容
多様体Ｍ上での最小値近傍ではＭに沿ってのエネルギー
関数が充分２次形式であると見なすことができる場合に
は、この近似は厳密の計算結果と一致するので、Ｓｍit
hのアルゴリズムと同等の正当性を持ち、かつ計算時間
は格段に速くなる。

【００８０】本実施形態のアルゴリズムによる問題の解
決方法について図１に従って述べる。本実施形態のアル
ゴリズムは以下のような計算手順を踏む。

【００８１】ステップS1で初期位置設定をx＝x₀とす
る。

【００８２】ステップS2で次の操作を行う。

【００８３】x:＝x+Appr_M(x) (40) ステップS3で、エネルギーの収束を判断する。以下、ス
テップS2〜S8を収束するまで行う。

【００８４】ステップS4でＧ_new:＝Proj_xGrad_xE (41) を計算する。

【００８５】ステップS5で反復回数を判定する。初回で
あれば、ステップS6に移行し、Ｇ:＝Ｇ_new、Ｈ:＝−Ｇ
_newとする。

【００８６】次に、ステップS8で１）曲線Ｃを

【外３９】

【００８７】とし、Ｃ上で、E(t)を最小とするtを求め
る。２）ベクトルGを測地線に沿ってtだけ近似的に平行移動
させたものを τ(Ｇ):＝Proj_x(Ｇ−tΓ(Ｈ，Ｇ)) (42) とする。３）ベクトルＨを測地線に沿ってtだけ近似的に平行移
動させたものを τ(Ｈ):＝Proj_x(Ｈ−tΓ(Ｈ，Ｈ)) (43) とする。４）点xを測地線上をtだけ移動させる。

【００８８】

【外４０】

【００８９】ステップS2に戻って、許容多様体Ｍへの移
動をさせた後に、ステップS3で収束判定を行う。判定
は、例えば、エネルギー差が所望の値より大きいか小さ
いかで判定する。充分収束したと判断した場合はプログ
ラム（サブルーチン、関数）を終了させる。

【００９０】ステップS5でＧ_newを計算した後、ステッ
プS6の判定を経て、ステップS8へ移行し、

【外４１】

【００９１】

【外４２】

【００９２】

【外４３】

【００９３】を計算代入してステップS8へ移行する。こ
の操作を収束するまで行う。

【００９４】次に、図１のフローチャートにつき説明し
た最適化方法を実現するための計算機の構成を説明をす
る。本実施形態の計算方法は、図５、６に示すように計
算機本体１１１とディスプレイ１１２、キーボード１１
３とマウス１１４とからによって構成される計算機シス
テム上とフロッピー（登録商標）ディスク１１０等の記
憶媒体によって実現されている。計算機本体内１１１に
は、CPU１１６、ハードディスク１１７、メモリ１１
８、フロッピーディスクドライブ１１３等が設置されて
いる。

【００９５】コンピュータは、メモリ１１８上にハード
ディスク１１７またはフロッピー・ディスク１１０等の
記憶媒体から本発明のアルゴリズムに従う最適化方法に
関する実行プログラムをロードし、必要な計算領域をメ
モリ１１８上に確保し、キーボード１１３から入力され
たり、フロッピーディスク１１０から読み込まれるなど
して適当な方法で入力されたエネルギー関数や拘束条件
を表わす関数と、初期位置に対して、図１のフローチャ
ートに対応する手順で計算を行う。得られた関数の最小
値や最小値を与える位置が数値あるいはグラフとしてデ
ィスプレイ１１２に表示される。

【００９６】〔実施例１〕実施例１としては以下の幾何
学的な問題を考察する。

【００９７】考える空間はn次元実ベクトル空間Ｒⁿ(e_i|
i=1，・・・，n)として、規格直交基底を(e_i|i=1，・・・，n)
とし、つまりb_ij=<e_i， e_j>=δ(i，j=1，・・・，n)とす
る。許容集合Ｍを

【外４４】

【００９８】また、エネルギー汎関数を

【外４５】

【００９９】とする。この場合、拘束条件を表わす条件
方程式としては

【外４６】

【０１００】を取ればよい。

【０１０１】便宜上のため、ここでeⁱ:=b^ije_jという双
対基底を導入すると

【外４７】

【０１０２】

【外４８】

【０１０３】となる。ヘシアンも次のように計算され
る。

【０１０４】 Hess_xS¹(U，V)＝2δ_ijUⁱV^j （50） Hess_xS²(U，V)＝0 これらの量を利用して、図１のアルゴリズムに従って計
算を行うと、図２のように収束する。図２の縦軸は、図
１のステップＳ３で評価したエネルギー関数Ｅ(x)の値
と、ステップＳ３からＳ８までの処理を行った後に評価
したエネルギー関数Ｅ(x)との値の差であり、横軸はそ
の回数である。但し、一度目は初期位置でのエネルギー
関数Ｅ(x)の値自身を与えた。図２より判るように、エ
ネルギーＥ(x)は収束した。

【０１０５】〔実施例２〕実施例２として、平面上の液
滴形状の問題に上述した実施形態の最適化アルゴリズム
を適応してみる。図３は、平面上の液滴を示す図であ
る。拘束条件として、液滴の体積は一定であり、表面積
を最小とするものとする。液滴３０１は図３に示すよう
に平面板３０２にの上に形成されている。液滴３０１に
回転対称を仮定し、更に半球形状と同じトポロジーを持
つことを仮定して、図３のような半径をパラメータとし
て、計算を行う。

【０１０６】拘束条件としては、液滴の体積を一定にす
るもので、その条件を満たした許容多様体は次の式で与
えられる。

【０１０７】Ｍ＝｛（r¹，r²，・・・，rⁿ）∈Ｒⁿ|S(r)=0｝ (51) ここで、

【外４９】

【０１０８】但し、

【外５０】

【０１０９】として、

【外５１】

【０１１０】としている。

【０１１１】また、エネルギー汎関数としては、

【外５２】

【０１１２】を考えている。

【０１１３】この時、上記アルゴリズムを使用するため
の幾何学量として、以下のような量が具体的に求まる。
但し、条件式が一つのため、S(r)をS¹(r)と同一視して
いる。

【０１１４】

【外５３】

【０１１５】

【外５４】

【０１１６】

【外５５】

【０１１７】

【外５６】

【０１１８】

【外５７】

【０１１９】

【外５８】

【０１２０】

【外５９】

【０１２１】例えば、S(r)のヘシアンは次のようにもと
まる。

【０１２２】 Hess_xS¹(U，V)＝2s_ijUⁱV^j （55）但し、

【外６０】

【０１２３】また、個々の項は以下で与えられている。

【０１２４】

【外６１】

【０１２５】

【外６２】

【０１２６】

【外６３】

【０１２７】としている。

【０１２８】同様に、エネルギー汎関数についても、エ
ネルギー汎関数を

【外６４】

【０１２９】と書き換えることにより、その方向微分等
々は以下のように計算される。

【０１３０】

【外６５】

【０１３１】

【外６６】

【０１３２】

【外６７】

【０１３３】

【外６８】

【０１３４】

【外６９】

【０１３５】同様にエネルギー汎関数のヘシアンHess_xE
(U，V)＝2e_ijUⁱV^jについても以下のように計算される。

【０１３６】

【外７０】

【０１３７】

【外７１】

【０１３８】

【外７２】

【０１３９】

【外７３】

【０１４０】

【外７４】

【０１４１】

【外７５】

【０１４２】これらの量を利用して、図１のアルゴリズ
ムに従って計算を行うと図４のように収束する。図４
は、適当な形状から出発して、計算を行った後に収束し
た際の液滴の形状を示したもので、厳密な解と一致して
いることがわかる。

【０１４３】なお、上述したように、本実施形態のアル
ゴリズムをサブルーチンとして、図６のフロッピーディ
スク１１０またはハードディスク１１７等のメディアに
収めておくこともできる。このことにより容易にシミュ
レーションができることとなる。但し、サブルーチンの
引数として、エネルギー関数や条件関数を取るとしてい
る。

【０１４４】上記実施形態では、実数体上のn次元ユー
クリッド空間において、拘束条件が独立なm個の等式で
指定される場合について、２次の精度で測地線方程式の
解を具体的に構築し、許容集合から離れた場合の補正を
考慮したアルゴリズムを提案した。

【０１４５】上記実施形態の方法を用いると、測地線方
程式を厳密に解く場合に比べて、計算量が著しく少なく
なる。一方、測地線方程式の厳密解に従って計算した場
合、現実には、数値誤差のために２次近似の場合と同
様、許容集合への回帰処理が必要になる。従って、実際
上、２次近似による計算は、厳密解を用いた計算と比較
して遜色がなく、計算量が少ない分、有利であると考え
られる。

【０１４６】また、上記実施形態で述べたアルゴリズム
においては、測地線の計算及びベクトルの平行移動の計
算において２次関数近似を行ったが、この内、どちらか
あるいは両方を１次近似に置き換えても、同様のアルゴ
リズムを構成でき、同様の計算結果を得ることができ
る。

【０１４７】なお、本発明は、上記実施形態の機能が実
現できる範囲において、複数の機器（例えばコンピュー
タ本体、インターフェイス機器、ディスプレイなど）か
ら構成されるシステムに適用しても、単一の機器からな
る装置に適用してもよい。

【０１４８】また、前述した実施形態の機能を実現する
ように各種デバイスを動作させることを目的として、該
各種デバイスと接続された装置あるいはシステム内のコ
ンピュータに、前述した実施形態の機能を実現するソフ
トウェアのプログラムコードを供給し、供給されたプロ
グラムにしたがって、そのシステムあるいは装置のコン
ピュータ（またはＣＰＵやＭＰＵ）により、前記各種デ
バイスを動作させることによって実施したものも、本願
発明の範囲に含まれる。またこの場合、記憶媒体から読
み出されたプログラムコード自体が前述した実施形態の
機能を実現することになり、そのプログラムコード自
体、そのプログラムコードをコンピュータに供給する手
段、例えばかかるプログラムコードを記憶した記憶媒体
は、本発明を構成することになる。

【０１４９】かかるプログラムコードを供給する為の記
憶媒体としては、例えば、フロッピーディスク、ハード
ディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯ
Ｍ、ＣＤ−Ｒ、磁気テープ、不揮発性のメモリカード、
ＲＯＭ等を使用できる。

【０１５０】また、コンピュータが読み出したプログラ
ムコードを実行することにより、前述した実施形態の機
能が実現されるだけでなく、そのプログラムコードの指
示に基づき、コンピュータ上で稼動しているＯＳ（オペ
レーティングシステム）、あるいは他のアプリケーショ
ンソフトなどと協働して前述の実施形態の機能が実現さ
れる場合にも、かかるプログラムコードは本願発明の範
囲に含まれることは言うまでもない。

【０１５１】更に、記憶媒体から読み出されたプログラ
ムコードが、コンピュータに挿入された機能拡張ボード
やコンピュータに接続された機能拡張ユニットに備わる
メモリに書き込まれた後、そのプログラムコードの指示
に基づき、その機能拡張ボードや機能拡張ユニットに備
わるＣＰＵ等が実際の処理の一部または全部を行い、そ
の処理によって前述した実施形態の機能が実現される場
合も含まれることは言うまでもない。

【０１５２】本願発明を上記記憶媒体に適用する場合、
その記憶媒体には、先に説明したフローチャートに対応
するプログラムコードを格納すればよい。

【０１５３】

【発明の効果】以上説明したように、本発明の最適化方
法を用いると簡単な計算手順で拘束条件付き最適化問題
の解を求めることができる。

【図面の簡単な説明】

【図１】拘束条件付き最適化問題の解を求める手順を示
すフローチャートである。

【図２】収束の様子を示す図である。

【図３】平面上の液滴を示す図である。

【図４】収束した際の液滴の形状を示した図である。

【図５】最適化方法を実行する計算機の外観構成を示す
図である。

【図６】最適化方法を実行する計算機のブロック構成を
示す図である。

【図７】従来の拘束条件付き最適化問題の解を求める手
順を示すフローチャートである。

Claims

【特許請求の範囲】

【請求項１】等式で与えられる拘束条件の下での最適
化関数の最小値を求める拘束条件付き最適化方法におい
て、条件式で定義される集合を有限次元実ベクトル空間内の
リーマン多様体と見なし、実ベクトル空間内の初期位置
から、拘束条件を満たす当該リーマン多様体へ接近する
接近工程と、前記リーマン多様体上の接ベクトルに対する測地線方程
式に関する指数関数写像を有限次数で打ち切り、該打ち
切りによって得られた近似測地線を実ベクトル空間内の
１次元軌道として生成する軌道生成工程と、前記リーマン多様体上の接ベクトルと前記軌道生成工程
によって生成される軌道に対して当該接ベクトルの平行
移動に関する指数関数写像の有限次数近似によって近似
平行移動を行う近似平行移動工程とを有することを特徴
とした拘束条件付き最適化方法。
【請求項２】前記接近工程によって前記リーマン多様
体に接近した地点において、最適化関数の勾配を前記実
ベクトル空間内の点として計算し、得られた勾配ベクト
ルを当該リーマン多様体の接空間に射影して接ベクトル
を更新生成する接ベクトル更新生成工程と、該接ベクトル更新生成工程により得られた接ベクトルに
基づいて、探索方向を表わす探索方向ベクトルを生成す
る探索方向ベクトル生成工程と、前記軌道生成工程によって生成された軌道に沿って、最
適化関数を最小化する軌道最小点を探索する軌道上最小
点探索工程とを備え、前記軌道生成工程では、前記探索方向ベクトルに基づい
て軌道が生成され、前記近似平行移動工程は、前記接ベクトルと、前記探索
方向ベクトルと、前記軌道最小点とに基づいて行なわ
れ、前記軌道最小点を前記接近工程において前記初期位置に
代えて用いることで前記各工程を繰り返すことを特徴と
する請求項１に記載の拘束条件付き最適化方法。
【請求項３】前記各工程を有限回繰り返し実行後に前
記接近工程によって前記リーマン多様体上に接近させた
点を、求める最小値を与える点とすることを特徴とする
請求項２に記載の拘束条件付き最適化方法。
【請求項４】前記接近工程によって前記リーマン多様
体上に接近させる毎に、当該接近させた点において前記
最適化関数の収束を判定する収束判定工程を備え、該収
束判定工程により収束したと判定された場合に、当該接
近させた点を求める最小値を与える点とすることを特徴
とする請求項３に記載の拘束条件付き最適化方法。
【請求項５】前記探索方向ベクトル生成工程では、前
記繰り返しの最初において、接ベクトル更新生成工程に
より得られた接ベクトルGに対して、前記探索方向ベク
トルHをH＝-Gにより生成することを特徴とする請求項２
に記載の拘束条件付き最適化方法。
【請求項６】前記探索方向ベクトル生成工程では、前
記繰り返しの２回目以降において、前回の接ベクトルを
G、前回の探索方向ベクトルをH、前記近似平行移動工程
により、接ベクトルGを近似平行移動したベクトルをτ
(G)、探索方向ベクトルHを近似平行移動したベクトルを
τ(H)、前記接ベクトル更新生成工程により更新された
接ベクトルをG'とするとき、新たな探索方向ベクトルH'
を【外１】により生成することを特徴とする請求項２に記載の拘束
条件付き最適化方法。
【請求項７】前記近似平行移動工程において、平行移
動に関する指数関数写像の有限次数近似が１次近似であ
ることを特徴とする請求項１に記載の拘束条件付き最適
化方法。
【請求項８】前記軌道生成工程において、測地線方程
式に関する指数関数写像の有限次数近似が２次近似であ
ることを特徴とする請求項１に記載の拘束条件付き最適
化方法。
【請求項９】コンピュータにより実行可能なプログラ
ムであって、請求項１に記載の最適化方法をコンピュー
タに実行させるためのプログラム。