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

Question

弾性体の変形を見るために，停留ポテンシャル法を用いようとしているのですが，その際に目的関数にひずみエネルギーをとり
独立な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;//探索方向再設定
}

echoes_x86 · Accepted Answer

おはようございます．

私の読解能力が無いのかも知れませんが…良く分からなくなってます．
論点を整理しますと，私は非線型関数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は定数として扱う)，
元の問題のものではありません．

色々と書きましたが，文献を読むなら，いきなり「参考にして」ではなく，
まずは「書いてあるままに」実装することをおすすめします．
何らかの改良を加えたい様子ですが，
「自分が何をしているか明解に説明できない」ものは「改良」とは言いません．
うまく動作しないならなおさらです．

「回答する」といいながら非常に申し訳ないのですが，
一度閉じることを検討した方がいいのではないかと思います．

echoes_x86 · Answer

こんにちは．

いくつかの誤解があるように思います．
>共役勾配法のみで
共役勾配法のみで計算を行うなら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+α)を通る下に凸の二次関数を計算しながら，
その最小値を予測するものです．
詳しく聞きたい場合は再び補足してください．

echoes_x86 · Answer

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

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

おはようございます．

こんにちは．

この回答への補足

こんばんは．

この回答への補足

似たような質問が見つかりました

関連するカテゴリからQ&Aを探す

デイリーランキングこのカテゴリの人気デイリーQ&Aランキング

マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング