sage_example -- Sage

Notebookの特徴

Notebookの特徴は、Notebookに説明文や数式を簡単に挿入できることです。

$$y = sin(2 \pi x) + \mathcal{N}(0,0.3)$$

# PRMLのsin曲線のデータ data = matrix([ [0.000000, 0.349486], [0.111111, 0.830839], [0.222222, 1.007332], [0.333333, 0.971507], [0.444444, 0.133066], [0.555556, 0.166823], [0.666667, -0.848307], [0.777778, -0.445686], [0.888889, -0.563567], [1.000000, 0.261502], ]); X = data.column(0) t = data.column(1) M = 3 # データのプロット x = var('x'); sin_plt = plot(sin(2*pi*x),[x, 0, 1], rgbcolor='green') data_plt = list_plot(zip(X, t)); data_plt (data_plt + sin_plt).show()

Sin曲線の多項式フィッティング

@interactコマンドでインタラクティブに変数の値を変更し、結果をみることができます。

# Φ関数定義 def _phi(x, j): return x^j

# PRMLの図1.3 @interact def _(M= slider(range(10))): # 計画行列Φ Phi = matrix([[ _phi(x,j) for j in range(0, (M+1))] for x in X.list()]) Phi_t = Phi.transpose() # ムーア＿ベンローズの疑似逆行列 Phi_dag = (Phi_t * Phi).inverse() * Phi_t; # 平均の重み Wml = Phi_dag * t f = lambda x : sum(Wml[i]*x^i for i in range(0, (M+1))) y_plt = plot(f, [x, 0, 1], rgbcolor='red') (y_plt + data_plt + sin_plt).show(ymin=-1.5, ymax=1.5)

Click to the left again to hide and once more to show the dynamic interactive window

共役勾配法

以下のような２次形式の関数を考えます。 $$ f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - 6 x_1 - 7 x_2 $$ 極値に達するには、勾配▽fからある程度接線方向tにずれた共役勾配$d_n$方向に進みます。 $$ d_n = - \nabla f(x_n) + \beta_n d_{n-1} $$

$\beta_n$は $$ \beta_n = \frac{(\nabla f(x_n))^T \nabla f(x_n)}{(\nabla f(x_{n-1}))^T \nabla f(x_{n-1})} $$ となり、dの初期値は$d_0 = - \nabla f(x_0)$から始めます。

$x$は刻み値$\alpha$、 $$ \alpha_n = - \frac{d_n^T \nabla f(x_n)}{d_n^T H d_n} $$ を使って次式で更新します。ここでHは、f(x)のヘッセ行列です。 $$ x_{n+1} = x_n + \alpha_n d_n $$

関数ｆの定義

# fを定義 def f(v): return 3/2 * v[0]^2 + v[0]*v[1] + v[1]^2 - 6*v[0] - 7*v[1]

変数のベクトル表現

# 変数定義 vars = var('x1 x2') v = vector([x1, x2])

▽ｆの計算

▽ｆの計算は、ちょっとトリックを使います。あらかじめ関数fの各変数での偏微分をdfsに保持しておき、その結果に引数のベクトル vxの値を代入した結果を返しています。

# 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])

ヘッセ行列

ヘッセ行列もSageの数式機能を使えば、簡単にもとめることができます。

# ヘッセ行列 H = matrix([[diff(diff(f(v),x_i), x_j) for x_i in v] for x_j in v]); print jsmath(H)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
3 & 1 \\
1 & 2
\end{array}\right)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr}
3 & 1 \\
1 & 2
\end{array}\right)

$\alpha_n$の定義

# α_nの定義 def alpha_n(x, d): return -d.dot_product(nabla_f(x)) / (d * H * d)

共役勾配法の反復処理

共役勾配法の反復処理は、至って単純です。条件を満たすまで与えられた式でxと共役勾配を更新するだけです。

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_n(x, d)*d if ((x - o_x).norm() < eps): break beta = nabla_f(x).dot_product(nabla_f(x)) / o_nabla_f_sqr d = -nabla_f(x) + beta*d if (d.norm() == 0): # 0割り対策 break k += 1 print "x=", x print "k=", k

x= (1, 3)
k= 2

x= (1, 3)
k= 2

Sageの最適化機能で結果を検証

求まった解x= (1, 3)をSageの最適化機能で求めた結果と比較します。

# 同様の処理をsageの機能を使って計算してみる g = 3/2*x1^2 + x1*x2 + x2^2 - 6*x1 -7*x2 minimize(g, [2, 1], algorithm="cg")

Optimization terminated successfully.
         Current function value: -13.500000
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
(1.0, 3.0)

Optimization terminated successfully.
         Current function value: -13.500000
         Iterations: 2
         Function evaluations: 5
         Gradient evaluations: 5
(1.0, 3.0)

# 関数と解をプロット p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4]) pt = point([1, 3, f(x)], color='red') (p3d+pt).show()

Sleeping...
If no image appears re-execute the cell. 3-D viewer has been updated.

sage_example

4947 days ago by takepwave