sage/PRML- 共役勾配法
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/22/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageで共役勾配法を試す [#ke371efb]
PRMLの第5章のニューラルネットワークの「勾配降下最適化」で...
方式がでてきたので、学習がてらsageを使って試してみます。
** 共役勾配法 [#hdaa5489]
[[朱鷺の杜Wiki>http://ibisforest.org/index.php]]の共役勾...
とプログラムの変数を合わせると、以下のような2次形式の関数...
$$
f(x) = \frac{1}{2} x^T A X - b^T x
$$
この時、極値に達するには、勾配▽fからある程度接線方向tにず...
方向に進まなくてはならない。（収束がほぼ直角に折れながら...
&ref(cg.png,,40%);
従って、
$$
d_n = - \nabla f(x_n) + \beta_n d_{n-1}
$$
と書けます（一つ前の\(d_{n-1}\)が接線方向tであることに注...
$$
\beta_n = \frac{(\nabla f(x_n))^T \nabla f(x_n)}{(\nabla ...
$$
となり、dの初期値は$d_0 = - \nabla f(x_0)$から始めます。
\(x\)は刻み値\(\alpha\)、
$$
\alpha_n = - \frac{d_n^T \nabla f(x_n)}{d_n^T A d_n}
$$
を使って次式で更新します。
$$
x_{n+1} = x_n + \alpha_n d_n
$$
sageへの入力:
#pre{{
# 共役勾配法(conjugate gradient method)
# f(x) = 1/2 x^T A x - b^T x
# d_0 = -▽f(x_0)
# d_n = -▽f(x_n)+(▽f(x_n)^T ▽f(x_n))/(▽f(x_n_1)^T ▽f(x_n_...
# x_n_p1 = x_n + t_n d_n
# t_n = - (d_n^T ▽f(x_n))/(d_n^R A d_n)
}}
** 関数、変数の定義 [#n4fb4edc]
syou6162さんの[[最適化理論][R]共役勾配法を実装してみた>ht...
の例題をSageを使って試してみます。
数式処理システムSageの特徴を活かすため、関数の引数をベク...
ベクトルvにセットします。
つぎに関数fを以下のように定義します。
$$
f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - t x_1 - 7 x_2
$$
sageへの入力:
#pre{{
# 変数定義
vars = var('x1 x2');
v = vector([x1, x2]);
}}
sageへの入力:
#pre{{
# 参考サイトhttp://d.hatena.ne.jp/syou6162/20090926/12539...
# 例題の関数：f = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2
# fを定義
def f(v):
return 3/2 * v[0]^2 + v[0]*v[1] + v[1]^2 - 6*v[0] - 7...
}}
*** ▽ｆの計算 [#z0bfddc9]
▽ｆの計算は、ちょっとトリックを使います。あらかじめ関数fの
各変数での偏微分をdfsに保持しておき、その結果に引数のベク...
vxの値を代入した結果を返しています。
sageへの入力:
#pre{{
# fを偏微分したリスト
dfs = [diff(f(v), x_i) for x_i in v];
# ▽fを定義 (dfsにvxの要素の値を適応した結果を返す)
def nabla_f(vx):
# ベクトルvxの各要素の値をvの要素に対応づける
s = dict(zip(v, vx));
# ベクトルの各要素の偏微分の結果にsを適応させる
return vector([df.subs(s) for df in dfs]);
}}
*** ヘッセ行列 [#p09285f2]
ヘッセ行列もSageの数式機能を使えば、簡単にもとめることが...
sageへの入力:
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for ...
print jsmath(H);
}}
sageの出力:
#pre{{
3 & 1 \\
1 & 2
}}
sageへの入力:
#pre{{
# α_kの定義
def alpha_k(x, d):
return -d.dot_product(nabla_f(x)) / (d * H * d);
}}
** 共役勾配法の反復処理 [#lbcf409f]
共役勾配法の反復処理は、至って単純です。条件を満たすまで...
更新するだけです。
sageへの入力:
#pre{{
# 共役勾配法の反復処理
eps = 0.001;
x0 = vector([2, 1]);
d = - nabla_f(vx=x0);
x = x0;
k = 1;
while (true):
o_nabla_f_sqr = nabla_f(x).dot_product(nabla_f(x));
o_x = x;
x += alpha_k(x, d)*d;
if ((x - o_x).norm() < eps):
break;
beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f...
d = -nabla_f(x) + beta*d;
if (d.norm() == 0): # 0割り対策
break;
k += 1;
print "x=", x;
print "k=", k;
}}
sageの出力:
#pre{{
x= (1, 3)
k= 2
}}
** Sageの最適化機能で結果を検証 [#x2f759cb]
求まった解x= (1, 3)をSageの最適化機能で求めた結果と比較し...
sageへの入力:
#pre{{
# 同様の処理をsageの機能を使って計算してみる
g = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2;
minimize(g, [2, 1], algorithm="cg")
}}
sageの出力:
#pre{{
Optimization terminated successfully.
Current function value: -13.500000
Iterations: 2
Function evaluations: 5
Gradient evaluations: 5
(1.0, 3.0)
}}
** 解のプロット [#d6a0206d]
Sageの強みは簡単に結果をグラフ化できることです。
sageへの入力:
#pre{{
# 関数と解をプロット
p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4]);
pt = point([1, 3, f(x)], color='red');
show(p3d+pt);
}}
&ref(1.png);
sageへの入力:
#pre{{
# 初期値を変えて収束の様子をプロットしながら再計算
pts = pt
pt2s = point([1, 3], color='red')
eps = 0.001
x0 = vector([-0.5, 1])
d = - nabla_f(vx=x0)
x = x0
k = 1
while (true):
o_nabla_f_sqr = nabla_f(x).dot_product(nabla_f(x))
o_x = x
pts += point([o_x[0], o_x[1], f(o_x)], color='yellow')
pt2s += point([o_x[0], o_x[1]], color='blue')
x += alpha_k(x, d)*d
if ((x - o_x).norm() < eps):
break
beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f...
d = -nabla_f(x) + beta*d
if (d.norm() == 0): # 0割り対策
break
k += 1
pts += point([x[0], x[1], f(x)], color='red')
pt2s += point([x[0], x[1]], color='red')
(p3d + pts).show()
}}
&ref(2.png);
** コンター図で解の収束を見る [#g4d38f68]
共役勾配法での解の収束の様子をコンターズで示します。
[[朱鷺の杜Wiki>http://ibisforest.org/index.php]]の共役勾...
の通り、2回の計算で収束しているのが分かります。
sageへの入力:
#pre{{
# コンターズで共役勾配と接線の法線方向の関係、および2回で...
cnt_plot = contour_plot(g, [x1, -1, 4], [x2, -1, 4], fill...
(cnt_plot + pt2s).show()
}}
&ref(sage0.png);
** コメント [#ofbfb6b1]
#vote(おもしろかった[5],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- pt2未定義のエラーを修正しました。修正箇所は、pt2s = poi...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/22/
また、Sageのサーバを公開しているサイト（http://sage.math....
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageで共役勾配法を試す [#ke371efb]
PRMLの第5章のニューラルネットワークの「勾配降下最適化」で...
方式がでてきたので、学習がてらsageを使って試してみます。
** 共役勾配法 [#hdaa5489]
[[朱鷺の杜Wiki>http://ibisforest.org/index.php]]の共役勾...
とプログラムの変数を合わせると、以下のような2次形式の関数...
$$
f(x) = \frac{1}{2} x^T A X - b^T x
$$
この時、極値に達するには、勾配▽fからある程度接線方向tにず...
方向に進まなくてはならない。（収束がほぼ直角に折れながら...
&ref(cg.png,,40%);
従って、
$$
d_n = - \nabla f(x_n) + \beta_n d_{n-1}
$$
と書けます（一つ前の\(d_{n-1}\)が接線方向tであることに注...
$$
\beta_n = \frac{(\nabla f(x_n))^T \nabla f(x_n)}{(\nabla ...
$$
となり、dの初期値は$d_0 = - \nabla f(x_0)$から始めます。
\(x\)は刻み値\(\alpha\)、
$$
\alpha_n = - \frac{d_n^T \nabla f(x_n)}{d_n^T A d_n}
$$
を使って次式で更新します。
$$
x_{n+1} = x_n + \alpha_n d_n
$$
sageへの入力:
#pre{{
# 共役勾配法(conjugate gradient method)
# f(x) = 1/2 x^T A x - b^T x
# d_0 = -▽f(x_0)
# d_n = -▽f(x_n)+(▽f(x_n)^T ▽f(x_n))/(▽f(x_n_1)^T ▽f(x_n_...
# x_n_p1 = x_n + t_n d_n
# t_n = - (d_n^T ▽f(x_n))/(d_n^R A d_n)
}}
** 関数、変数の定義 [#n4fb4edc]
syou6162さんの[[最適化理論][R]共役勾配法を実装してみた>ht...
の例題をSageを使って試してみます。
数式処理システムSageの特徴を活かすため、関数の引数をベク...
ベクトルvにセットします。
つぎに関数fを以下のように定義します。
$$
f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - t x_1 - 7 x_2
$$
sageへの入力:
#pre{{
# 変数定義
vars = var('x1 x2');
v = vector([x1, x2]);
}}
sageへの入力:
#pre{{
# 参考サイトhttp://d.hatena.ne.jp/syou6162/20090926/12539...
# 例題の関数：f = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2
# fを定義
def f(v):
return 3/2 * v[0]^2 + v[0]*v[1] + v[1]^2 - 6*v[0] - 7...
}}
*** ▽ｆの計算 [#z0bfddc9]
▽ｆの計算は、ちょっとトリックを使います。あらかじめ関数fの
各変数での偏微分をdfsに保持しておき、その結果に引数のベク...
vxの値を代入した結果を返しています。
sageへの入力:
#pre{{
# fを偏微分したリスト
dfs = [diff(f(v), x_i) for x_i in v];
# ▽fを定義 (dfsにvxの要素の値を適応した結果を返す)
def nabla_f(vx):
# ベクトルvxの各要素の値をvの要素に対応づける
s = dict(zip(v, vx));
# ベクトルの各要素の偏微分の結果にsを適応させる
return vector([df.subs(s) for df in dfs]);
}}
*** ヘッセ行列 [#p09285f2]
ヘッセ行列もSageの数式機能を使えば、簡単にもとめることが...
sageへの入力:
#pre{{
# ヘッセ行列
H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for ...
print jsmath(H);
}}
sageの出力:
#pre{{
3 & 1 \\
1 & 2
}}
sageへの入力:
#pre{{
# α_kの定義
def alpha_k(x, d):
return -d.dot_product(nabla_f(x)) / (d * H * d);
}}
** 共役勾配法の反復処理 [#lbcf409f]
共役勾配法の反復処理は、至って単純です。条件を満たすまで...
更新するだけです。
sageへの入力:
#pre{{
# 共役勾配法の反復処理
eps = 0.001;
x0 = vector([2, 1]);
d = - nabla_f(vx=x0);
x = x0;
k = 1;
while (true):
o_nabla_f_sqr = nabla_f(x).dot_product(nabla_f(x));
o_x = x;
x += alpha_k(x, d)*d;
if ((x - o_x).norm() < eps):
break;
beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f...
d = -nabla_f(x) + beta*d;
if (d.norm() == 0): # 0割り対策
break;
k += 1;
print "x=", x;
print "k=", k;
}}
sageの出力:
#pre{{
x= (1, 3)
k= 2
}}
** Sageの最適化機能で結果を検証 [#x2f759cb]
求まった解x= (1, 3)をSageの最適化機能で求めた結果と比較し...
sageへの入力:
#pre{{
# 同様の処理をsageの機能を使って計算してみる
g = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2;
minimize(g, [2, 1], algorithm="cg")
}}
sageの出力:
#pre{{
Optimization terminated successfully.
Current function value: -13.500000
Iterations: 2
Function evaluations: 5
Gradient evaluations: 5
(1.0, 3.0)
}}
** 解のプロット [#d6a0206d]
Sageの強みは簡単に結果をグラフ化できることです。
sageへの入力:
#pre{{
# 関数と解をプロット
p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4]);
pt = point([1, 3, f(x)], color='red');
show(p3d+pt);
}}
&ref(1.png);
sageへの入力:
#pre{{
# 初期値を変えて収束の様子をプロットしながら再計算
pts = pt
pt2s = point([1, 3], color='red')
eps = 0.001
x0 = vector([-0.5, 1])
d = - nabla_f(vx=x0)
x = x0
k = 1
while (true):
o_nabla_f_sqr = nabla_f(x).dot_product(nabla_f(x))
o_x = x
pts += point([o_x[0], o_x[1], f(o_x)], color='yellow')
pt2s += point([o_x[0], o_x[1]], color='blue')
x += alpha_k(x, d)*d
if ((x - o_x).norm() < eps):
break
beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f...
d = -nabla_f(x) + beta*d
if (d.norm() == 0): # 0割り対策
break
k += 1
pts += point([x[0], x[1], f(x)], color='red')
pt2s += point([x[0], x[1]], color='red')
(p3d + pts).show()
}}
&ref(2.png);
** コンター図で解の収束を見る [#g4d38f68]
共役勾配法での解の収束の様子をコンターズで示します。
[[朱鷺の杜Wiki>http://ibisforest.org/index.php]]の共役勾...
の通り、2回の計算で収束しているのが分かります。
sageへの入力:
#pre{{
# コンターズで共役勾配と接線の法線方向の関係、および2回で...
cnt_plot = contour_plot(g, [x1, -1, 4], [x2, -1, 4], fill...
(cnt_plot + pt2s).show()
}}
&ref(sage0.png);
** コメント [#ofbfb6b1]
#vote(おもしろかった[5],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- pt2未定義のエラーを修正しました。修正箇所は、pt2s = poi...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
