共役勾配法を用いた最適化について

質問

弾性体の変形を見るために，停留ポテンシャル法を用いようとしているのですが，その際に目的関数にひずみエネルギーをとり
独立な24変数の関数の極小値を共役勾配法で求めようとしています．

Hesse行列を手計算でだし，最適化を行うプログラムを書いたのですが
どうも上手くいきません．

下にプログラムの一部分をのせます．
共役勾配法について勘違いしている可能性もあるので
問題点を教えて頂ければ幸いです．

x_vectorは24変数ベクトルで，dは探索方向，HはHesse行列です．
for(i=0;i<10;i++)
{
lambda = (x_vector*d)/(d*H*d);

x_save_vector = x_vector;
x_vector=x_vector - lambda*H*d;//変数べクトル再設定

if(i==0)
{
// 初期値
gamma = 1;
}
else
{
gamma =x_vector*x_vector/(x_save_vector*x_save_vector);
}

d = x_vector + gamma*d;//探索方向再設定
}

回答No.3

echoes_x86 · ベストアンサー

おはようございます．

私の読解能力が無いのかも知れませんが…良く分からなくなってます．
論点を整理しますと，私は非線型関数F(x)の最小化に共役勾配法を直接適用するつもりで回答しています．
これは，質問者様の最初の補足がそうなっているからです．
> Newton法は使用せず、共役勾配法のみで解こうとしています。

Hessianを使わない「普通の」共役勾配法では，
F(x)を直接相手にするわけですから，
ある探索方向dが決まったとして，その方向に沿った切断面が二次関数になるとは限りません．
また，dに沿った切断面の評価関数の微分を計算したり(最小値候補が解析的に計算できます)，
値を調べるのに計算コストがかかる場合があります
(そもそも共役勾配法はこうした大規模問題に適用するものです)．
ですから，最低3点の値を評価して評価関数の切断面を二次関数であるものとし，
その最小値の予測を反復する必要が生じるわけです．
参考文献[2]にも「評価関数(f(x))が二次関数の時」，「厳密直線探索」すれば… と書いてあります．
これは二次の評価関数 x^T*A*x の場合に，
ある探索方向dに関する評価関数 (x+αd)^T*A*(x+αd) をαで微分して零として最小値を与えるαを解析的に計算すれば，
厳密直線探索ができるという流れです．
非線型関数のHessianを使うとはどこにも書いてないですよね？

ですが，質問者様は
・Hessianを利用する
・αは解析的に決める
ことにこだわっています．
ここから推測される状況はNewton法の更新方程式 H*d = -g と等価な
二次の最小化問題 d^T*H*d/2 + g'*d を解こうとしているとしか思えません(gは勾配)．
この場合はHessianが正定値になりさえすれば，αを解析的に求めて良いでしょう．
「二次近似した」評価関数の切断面が二次関数になると分かっているからです(つまり微分して零の点が唯一の最小値)．
Newton法のHessianが極端に疎だったり大規模だったりする場合に，
H*d = -gの求解に共役勾配法のような反復解法を用いることは考えられます．
この場合，Newton法の内部で利用される共役勾配法が収束して初めて元の問題の解を更新することになります．
共役勾配法の評価関数はd^T*H*d/2 + g'*dであり(ここでHとgは定数として扱う)，
元の問題のものではありません．

色々と書きましたが，文献を読むなら，いきなり「参考にして」ではなく，
まずは「書いてあるままに」実装することをおすすめします．
何らかの改良を加えたい様子ですが，
「自分が何をしているか明解に説明できない」ものは「改良」とは言いません．
うまく動作しないならなおさらです．

「回答する」といいながら非常に申し訳ないのですが，
一度閉じることを検討した方がいいのではないかと思います．

お礼

おはようございます。

＞非線型関数のHessianを使うとはどこにも書いてないですよね？
確かに正定値対象行列Aと書いてあるけど、Hesse行列とは
書いてないですね…。

私自身の理解が足りないと思うので、
参考書をゆっくり読んで、もう一回考え直してきます。

何回も丁寧な回答ありがとうございました。

回答No.2

echoes_x86

こんにちは．

いくつかの誤解があるように思います．
>共役勾配法のみで
共役勾配法のみで計算を行うなら24回云々は忘れてください．
その代わり勾配ベクタを反復ごとに再計算してください．
また，ステップ幅αと共役方向生成βに混同(？)が見られます．
質問者様のβはおそらくFletcher-Reeves公式であると思われますが，
αは私は見たことがありません．
無論私が見たことが無いから間違いだなどとは言いませんが，
このαは普通は直線探索で決定するものです．
一方，質問者様のαはdkに対してH-共役な方向を生成するように読めます．
再度指摘しますが，普通の共役勾配法ではHessianは用いません．
用いるものがあってもそれはβの計算で，それも毎回は更新しないというものだったと思います．
また，定期的な再出発(リスタート)が必要です．
これは何回かごとにβ=0として探索方向を勾配方向にとりなおすものです．

>目的関数の値が増大する
直線探索をする以外に解決案はないと思われます．
共役勾配法は「勾配を基準としつつも前の数回とは違う方向を向く(共役方向と言います)」
ことを指針とします．
ここで，「どちらに進むか」を決めるのがβ，「どれだけ進むのか」を決めるのがαです．
そして勾配(=局所的に最も減る方向)とはβだけ違う方向を向くというのが微妙で，
要するに目的関数が減少することを保証できません．
なので，その方向の最小値を計算するために，
「どれだけ進むのか」を毎回計算で求めます．
これがαに関する一次元最小化問題 min F(xk+α) です(直線探索)．
これは種々の方法がありますが，簡単にやるなら二次内挿でいいでしょう．
大雑把に言って以下のものです．
(1)F(xk) > F(xk+0.5α) かつ F(xk+0.5α) < F(xk+α)のαを見つける
(2)上記の範囲の最小値を与えるαを確定
これはF(xk)，F(xk+0.5α) ，F(xk+α)を通る下に凸の二次関数を計算しながら，
その最小値を予測するものです．
詳しく聞きたい場合は再び補足してください．

補足

こんばんわ。
丁寧な回答ありがとうございます。
βの決め方はFletcher-Reevesの公式から決めています。
αですが、私が参考にしている本[1]の共役方向法の所に書いてあったので
それを参考にきめました。
ネット上では[2]の共役方向法の部分と同じように決めました。
ひょっとしたら私は共役方向法を使っているのかもしれません。
最適化は独学なので勘違いが多い可能性は高いです。

直線探索について
二次内挿というのは２次補間法のことですよね？

今プログラムを直し直しシミュレーションをしているのですが
一応最初よりましになってきました。

[1]非線形計画法　著者:今野浩，山下浩
[2]http://www.az.cs.is.nagoya-u.ac.jp/jsiam/tutorial-shape-opt-4.pdf

回答No.1

echoes_x86

こんばんは．

コードはあまり読んでませんが…
気になる点があるので補足をお願いします．
まず何よりも「どうも上手くいきません」の内容を補足してください．
また，Hessianまで考えているということは非線型関数の最適化という意味での共役勾配法と思われますが，
・探索方向に勾配(最急降下方向)を計算する部分が無い
・直線探索関数が無い
の2点が気になります．
それとも，最適化の本命はNewton法で，
その更新方程式 H*d = -g を解くために線型共役勾配法を適用しよう，
ということでしょうか．
普通の非線型共役勾配法はHessianを陽に計算しないからです(一部例外もありますが)．
とりあえず，後者を仮定します．
この場合，10回以内に H*d = -g の解が収束することを期待していることになりますが，
変数の次元を考えるとこれは無理である気がします．
Hessianが正定値だとして，厳密な直線探索を行えば 24回以内の探索で収束が保証されますが，
Hessianの固有値の分布によっては非常に収束が遅くなりますし，直線探索無しではかなり厳しいでしょう．

補足

回答ありがとうございます。

「どうも上手くいきません」というのは、目的関数の値
が増大していくことです．

やろうとしていることは非線形関数の最適化です。

・探索方向に勾配(最急降下方向)を計算する部分が無い
・直線探索関数が無い
の二点ですが、私の不勉強でした。

今、参考書をよんでいるのですが、上の二つについて
確かに書いてありました。
さらに直線探索の回数についても記述してありました。
だからソースコードは間違ってます。すいません…

Newton法は使用せず、共役勾配法のみで解こうとしています。

直線探索関数については、まだあまりわかっていません。
とりあえず次のように探索していこうと考えています。

Δf(x_k):勾配ベクトル,　Q:Hesse行列

ステップ幅:α_k=-{Δf(x_k).d_k}/{(d_k)^T.Q.d_k}

探索方向での探索：x_k+1 = x_k + α_k　.d_k

探索を24回以上行い，探索方向を再決定

探索方向は次のように求めようとしています。

β_k+1={Δf(x_k+1).Δf(x_k+1)}/{Δf(x_k).Δf(x_k)}

d_k+1 = -Δf(x_k+1)　+　β_k+1.d_k