![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
sageは、mathematicaのような数式処理を行うオープンソースのソフトウェアです。sageの歴史はまだ浅く、ウィリアム・スタイン 氏(William Stein)によって2005年2月に開発がスタートし、2006年2月のUCSD SAGE Days 1でSage 1.0が公開され、最新のバージョンは4.7.2です（注１）。
Sageの特徴を挙げると、
があります。
Sageの最大の特徴は、 FirefoxやInternet Explorer等のブラウザーからSage Notebook Serverにアクセスして、気軽に数式処理を実行することが出来ることです（注２）。
Notebookは、Sageでの一連の計算を記録したノートであり、計算に関する説明文を挿入したり、値を変更して再計算することができます。
Sageのノートブックを体験するには、Sageの開発サイトでアカウントを作成し、ワークシートを作成するのが最も簡単な方法です（注３）。
初期画面には、1個のセルが表示され、ここにpython形式で数式を入力し、evaluateボタンで結果が表示されます。
ノートブックには、説明文や図、数式を入力するためのツールが提供されています。セルの上下にマウスを移動すると青い帯が表示されますので、シフトキーとクリックを同時に押すと図１のようなエディタ画面が表示されます。
説明文には、Latex形式で数式が書けるので、計算手法の解説などが簡単に挿入でき、ノートブックの活用に役立ちます。
&re(editor_screen.png);
図１ エディタ画面の編集例
Sageのとても便利な機能にグラフの重ね書きがあります（注４）。
例として、Sin曲線のフィッティング用データと元のSin曲線を重ね書きする処理はリスト１のようになります。
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曲線のフィッティング用データと元のSin曲線の重ね書き例
様々な値に対する結果をインタラクティブに提供するために、@interactコマンドが用意されています。
@interactに続いて、関数定義def _(引数)を記述します。
からユーザが値を指定することができます。
図２にSin曲線のフィッティング問題で多項式の次数Mに０から９までの値のリストを与えた例を示します。
図２ Sin曲線のフィッティング問題で多項式の次数Mをスライダーバーで与えた例
共役勾配法を例にSageの数式処理システムらしい解法を紹介します。 Sageの数式処理機能とPythonの記述力を合わせると、とてもスマートに共役勾配法で2次関数の極値を求めることができます。
以下のような２次形式の関数を考えます。
極値に達するには、勾配▽fからある程度接線方向tにずれた共役勾配dn方向に進みます。
βnは
となり、dの初期値はd0=−∇f(x0) から始めます。 xは刻み値α、
を使って次式で更新します。 ここでHは、f(x)のヘッセ行列です。
ベクトルvを変数x1, x2で定義し、これを使って関数fを定義します。
# 変数定義
vars = var('x1 x2')
v = vector([x1, x2])
# fを定義
def f(v):
return 3/2 * v[0]^2 + v[0]*v[1] + v[1]^2 - 6*v[0] - 7*v[1]
次に、▽ｆを計算します。 あらかじめ関数ｆの各変数での偏微分した結果を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)
αnの定義も式の通りです。#pre{{
# α_nの定義 def alpha_n(x, d):
return -d.dot_product(nabla_f(x)) / (d * H * d)
}}
すべての準備が整ったので、共益勾配法を使って極値を計算してみます。
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
となり、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)
求まった解をプロットすると、図３のようになります。
# 関数と解をプロット p3d = plot3d(g, [x1, -1, 4], [x2, -1, 4]) pt = point([1, 3, f(x)], color='red') (p3d+pt).show()
図３ 関数fと求めた解をプロットした図
Sageではすべての処理をSage単独で行うのではなく、既存のツールと連携するためのインタフェースを提供しています。
例として、Rでの主成分分析の結果をSageに渡す方法を紹介します（注５）。
Oil Flowのデータを使って主成分分析をします。データは、あらかじめノートブックにアップロードしておきます。 Rとのインタフェース関数rを使ってRのコマンドを記述します。
データを読み込み、主成分分析をresult変数に代入します。
# Rでデータを読み込みPCAを計算
fileName = DATA + 'DataTrn.txt'
oilflow = r("oilflow <- read.table('%s')" %fileName) result = r("result <- prcomp(oilflow)")
RのグラフをSageで表示するには、一度ファイルに出力する必要があります（図４参照）。
# ラベルの読み込み
fileName = DATA + 'DataTrnLbls.txt'
labels = r("oilflow.labels <- read.table('%s')" %fileName)
# プロットファイル名の設定
filename = DATA+'pca.pdf'
r.pdf(file='"%s"' %filename)
# 結果のプロット
r("col <- colSums(t(oilflow.labels) * c(4,3,2))")
r("pch <- colSums(t(oilflow.labels) * c(3,1,4))")
r("plot(result$x[,1:2], col=col, pch=pch, xlim=c(-3,3), ylim=c(-3,3))")
r.dev_off()
# 式を変えたときにはブラウザーで再読込必要 html('<img src="pca.pdf">')
図４ 主成分分析の結果をSageで表示した図
Rのデータをsageに渡しすためにsageobj関数と_save_メソッドが提供されています。
例として、Rで読み込んだデータセットをSageの3次元プロットを使ってプロットしてみます（図４参照）。
#Rのデータセットをsageの形式に変換すると'DATA'ディクショナリに列単位でV1, V2のようにセットされる lb = sageobj(labels)
# これをzipでまとめて使う lbs = zip(lb['DATA']['V1'],lb['DATA']['V2'],lb['DATA']['V3'])
# point3dでプロット
N = len(lbs)
plt = Graphics()
for n in range(N):
[x, y, z] = rs[n]
if lbs[n][0] == 1:
plt += point3d([x, y, z], rgbcolor='blue')
elif lbs[n][1] == 1:
plt += point3d([x, y, z], rgbcolor='green')
else:
plt += point3d([x, y, z], rgbcolor='red') plt.show()
図５ Rで読み込んだデータをSageでプロットした図
このように柔軟で便利な機能をもつSageですが、日本国内ではまだまだ普及していません。本稿を通じて、一人でも多くの方にSageを使って頂ければ幸いです。今後のユーザ会などの支援団体の創設や解説書の出版が待たれます。 沼田泰英氏がSageのレファレンスカードを日本語に翻訳して公開されています（注６）。印刷して傍らに置いておくと便利です。
