7. 勾配法のプログラミング(2)

7. 勾配法のプログラミング (2)
田中雅博
最適化プログラミング、教科書 pp.79–85
多変数関数の極値を数値的に求める代表的な手法その１。多変数の場合。
1
多変数の関数
2 変数関数 f = f (x, y) の最大値を求める問題を考える。1 変数の場合同様、
∂f
∂f
= 0,
=0
∂x
∂y
となるような (x, y) を求める問題であるが、解析的に解けないものとしよう。そのため、勾配法で
数値的に求める。
教科書 p.82 の図 3.2 を見よ。k 回目の点 (xk , yk ) から次の (xi+1 , yk+1 ) に進むところを考えよう。
(xk , yk ) で最も勾配の大きな方向は、この点での勾配ベクトル
(
)
∂f
∂x
∇fk =
∂f
∂y
x=xk ,y=yk
である。
そこで、この方向で関数が極大になる点を探す。これを、直線探索という。探索の幅の作り方は、
1 次元の場合と同じでよい（黄金探索法など、別の探索法もあるが、ここでは使用しない）
。
(
)
x
今後の便宜のために、x =
と置こう。なお、以下の議論では、x は 2 次元だけでなく、
y
一般に n 次元としてよい。
直線探索は、
(∇fk はベクトルである）
g(t) = f (xk + t∇fk )
とおく。2 次元ベクトルの場合、
g(t) = f (xk + t
∂f
∂f
, yk + t
)
∂xk
∂yk
となり、最大化は
max g(t)
t
と書ける。
procedure search(f(x),f’(x))
1. x の初期値を与え、h ← h0 とする。k = 0。
1
2. g(t) = f (xk + t∇fk ) に対して 1 次元探索を行い、g(t) を最大化する t を求める。
3. 得られた t を用いて、
∆xk = t∇fk , xk+1 = xk + ∆xk
とし、k = k+1。Step 2. からここまで、∥∆x∥ < δ となるまで繰り返す。
□
以下、一次元探索のアルゴリズムを再度示す。前回の 1 次元の問題のアルゴリズムと同じである
が、g は t の関数であることに注意する。
procedure linsearch(f(x),f’(x))
1. t の初期値を与え、h ← h0 とする。
2. 次のように置く。
h ← sgn(g ′ (t))|h|, T ← t,
T′ ← t + h
3. もし g(t) < g(T ′ ) であれば、次の計算を行う。
(a) g(t) ≥ g(T ′ ) となるまで次の計算を繰り返す。
h ← 2h,
T ← T ′,
T′ ← T + h
(b) t ← T, h ← h/2 と置く。
4. そうでなければ、次の計算を行う。
(a) g(t) ≤ g(T ′ ) となるまで、次の計算を繰り返す。
h←
h
,
2
T′ ← T′ − h
(b) t ← T ′ , h ← 2h と置く。
5. ステップ２に戻り、これを |g ′ (t)| ≤ ϵ となるまで繰り返す。
6. 得られた t を返す。
ここで、
∑ ∂f ∂fk
dg ∑ ∂f ∂xi
=
=
dt
∂xi dt
∂xi ∂xi
i
i
例題
f (x, y) = −(x − y)2 − (x − 1)2
の最大値を求める。
g(t) = −
((
) (
))2 ((
)
)2
∂f
∂f
∂f
xk + t
− yk + t
−
xk + t
−1
∂xk
∂yk
∂xk
2
∑ ∂f ∂fk
dg ∑ ∂f ∂xi
=
=
dt
∂xi dt
∂xi ∂xi
i
i
= −2 ((x − y) + (x − 1))
∂f
∂f
+ 2(x − y)
∂xk
∂yk
これの MATLAB プログラムを示そう。関数 g(t) を g( t,x_k,y_k,fx_k,fy_k )、その導関数
g ′ (t) を gp(t,x_k,y_k,fx_k,fy_k) とする。main 関数を hill_climbing(x_k,y_k)、線形探索
を linsearch(x_k,y_k) とする。
function hill_climbing(x_k,y_k)
delta2=1.0e-6;
delx=1;
dely=1;
count=1;
fx_k=-2*(x_k-y_k)-2*(x_k-1);
fy_k = 2*(x_k-y_k);
xs(1)=x_k;
ys(1)=y_k;
fx_k=-2*(x_k-y_k)-2*(x_k-1);
fy_k = 2*(x_k-y_k);
while (delx^2 + dely^2) > delta2
count=count+1
t = linsearch(x_k,y_k);
delx = t * fx_k;
dely = t * fy_k;
x_k = x_k + delx;
y_k = y_k + dely;
fx_k = -2 * (x_k - y_k) - 2 * (x_k -1);
fy_k = 2 * (x_k - y_k);
disp([x_k y_k])
xs(count)=x_k;
ys(count)=y_k;
end
’optimal solution’
disp([x_k, y_k])
plot(xs,ys,’o-’)
hold on
x=-2:0.2:2;
y=-2:0.2:2;
[xx,yy]=meshgrid(x,y);
z=kansu2(xx,yy);
contour(xx,yy,z)
3
end
以下、線形探索のプログラムである。
function t = linsearch(x_k,y_k)
t=0;
fx_k=-2*(x_k-y_k)-2*(x_k-1);
fy_k = 2*(x_k-y_k);
%pause
if gp(0,x_k,y_k,fx_k,fy_k)>0
h0=0.1;
else
h0=-0.1;
end
%h0
h=h0;
epsilon=1.0e-6;
while abs(gp(t,x_k,y_k,fx_k,fy_k)) > epsilon
h=sign(gp(t,x_k,y_k,fx_k,fy_k))*abs(h);
T = t;
Tp = t + h;
if g(T,x_k,y_k,fx_k,fy_k)<g(Tp,x_k,y_k,fx_k,fy_k)
’case1’
%pause
while g(T,x_k,y_k,fx_k,fy_k) < g(Tp,x_k,y_k,fx_k,fy_k)
h=h*2;
T=Tp;
Tp=T+h;
end
t=T;
h=h/2;
else
’case2’
while g(T,x_k,y_k,fx_k,fy_k) > g(Tp,x_k,y_k,fx_k,fy_k)
h=h/2;
Tp=Tp-h;
%pause
end
t=Tp;
h=2*h;
4
end
t
end
以下、関数 g(t) のプログラム。
function f = g( t,x_k,y_k,fx_k,fy_k )
x = x_k + t * fx_k;
y = y_k + t * fy_k;
f=-(x-y)^2-(x-1)^2;
end
以下は、g’(t) のプログラム。
function dg = gp(t,x_k,y_k,fx_k,fy_k)
x = x_k + t * fx_k;
y = y_k + t * fy_k;
dg = -2*((x-y)+(x-1))*fx_k + 2*(x-y)*fy_k
end
以下、もとの関数。
function z = kansu2( x,y )
z=-(x-y).^2 - (x-1).^2;
end
実際に試してみよう。
>> hill_climbing(-0.5,0)
0.3293
-0.2073
0.5573
0.7049
0.8021
0.6437
0.8694
0.9129
0.9416
0.8948
0.9614
0.9743
0.9828
0.9690
0.9886
0.9924
5
0.9949
0.9908
0.9966
0.9978
0.9985
0.9973
0.9991
0.9994
0.9996
0.9993
ans =
optimal solution
0.9996
0.9993
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2
-1.5
-1
-0.5
0
0.5
1
1.5
6
2
11 月 12 日課題
1. 上記で与えられた関数の最大化について、いくつかの点からスタートして、収束するかどう
か、試してみよ。
2. 次に、f (x, y) = −(cos x − x)2 − (x − y)4 の最大化を行ってみよ。収束状況について、言及
せよ。
[提出物] Word に、以下の内容を示せ。
• 学籍番号、氏名、課題出題日
• 問１について
– MATLAB プログラムリスト
– 収束した結果を示し、それについて考察せよ。
• 問２について
– MATLAB プログラムリスト
– プログラムの中で、問１とは違う部分をわかりやすく示せ（プログラムの中の該当す
る部分をマーカーで塗るか、四角で囲むなど）
– プログラムをそのように変更した理由、収束結果についての考察
• 感想
7

Download Report