![[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 初期ログイン時のノートブック画面
ノートブック画面でNew Worksheetをクリックすると新しいワークシートが作成されます。
ワークシートで式を評価するには、セルと呼ばれるテキストエリアを利用します。
セルの基本操作は、以下のように行います。
Sageの関数の使い方は、関数名のあとに?を付けてshift-returnでヘルプ画面が表示されます（図２参照）。
abs?
図２ abs関数のヘルプ画面
関数名やPythonの変数名は、タブキーで補完することができます。
例として、facを入力した後にタブキーを入力すると自動的にfactorまで入力が補完され、factorとfactorialの２つの候補が表示されています（図３）。
図３ facに対する補完結果
ワークシートには、説明文や図、数式を入力するためのツールが提供されています。セルの上下にマウスを移動すると青い帯が表示されますので、シフトキーとクリックを同時に押すと図４のようなエディタ画面が表示されます。
説明文には、Latex形式で数式が書けるので、計算手法の解説などが簡単に挿入でき、ワークシートの活用に役立ちます。
図４ エディタ画面の編集例
Sageでは、Pythonのインタプリタを使用しているため、Pythonの文法で記述します。Pythonでは#以下がコメント文となり、セミコロン;で複数の式を1行にまとめて記述することができます。
Sage固有の使い方としては、変数への代入結果を表示するためにセミコロンの後に変数名を追加します（図５の第1セル）。 Sageは数式処理システムですので、結果は数値ではなく数式を評価した値となります（図５の第2セル）。 数値として表示するには、N関数を使用します（図５の第3セル）
図５ Sage固有の使い方
Sageの数式で変数として認識させるには、変数をvar関数で宣言しなくてはなりません。変数の宣言はvarの引数に文字列で使用する変数名をスペースで区切って行います（図６の第1セル）。
Sageで関数を定義するには、以下の3つの方法があります。
例として3乗根をシンボリック関数、Python関数、lamba関数で定義し、その結果を表示してみます（図６の第3セル以降）。
図６ Sage変数と関数の定義例
Sageでよく使われる定数を表１に示します。
表１ Sageでよく使われる定数
| 円周率 | \( \pi \) | pi |
| 自然対数の底 | \( e \) | e |
| 虚数単位 | \( I \) | I |
| 無限大 | \( \infty\) | oo |
Sageのとても便利な機能にグラフの重ね書きがあります。
例として、Sin曲線＋ノイズのデータと元のSin曲線を重ね書きする処理を図７（第1セル）に示します。
図７ Sin曲線＋ノイズのデータ（100個）と元のSin曲線の重ね書き例
様々な値に対する結果をインタラクティブに提供するために、@interactコマンドが用意されています。
@interactに続いて、関数定義def _(引数)を記述します。
からユーザが値を指定することができます。
図８にSin曲線のフィッティング問題で多項式の次数Mに０から９までの値のリストを与えた例を示します。
図８ Sin曲線のフィッティング問題で多項式の次数Mをスライダーバーで与えた例
本稿で紹介したワークシートや参考文献[1]で使用したワークシート（注4）は、自由にダウンロードできます。
各公開ワークシートを表示し、先頭にあるDownloadを押すと、ダウンロードが開始します。 ダウンロードしたワークシートは、ユーザのノートブックでUploadをクリックし、 ダウンロードしたファイルを指定することで、簡単にアップロードできます。 自分のノートブックでワークシートの計算を再現したり、値を変えて評価することでができますので、是非試してみて下さい。
共役勾配法を例にSageの数式処理システムらしい解法を紹介します。 Sageの数式処理機能とPythonの記述力を合わせると、とてもスマートに共役勾配法で2次関数の極値を求めることができます。
以下のような２次形式の関数を考えます。 $$ f(x) = \frac{3}{2} x_1^2 + x_1 x_2 + x_2^2 - 6 x_1 - 7 x_2 （式１） $$
極値に達するには、勾配▽fからある程度接線方向tにずれた共役勾配dn方向に進みます。
$$ d_n = - \nabla f(x_n) + \beta_n d_{n-1} （式２） $$
β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の初期値はd0=−∇f(x0) から始めます。 xは刻み値α、 $$ \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 （式５） $$
ベクトル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)
\( \alpha_n \) の定義も式の通りです。
# α_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のレファレンスカードを日本語に翻訳して公開されています（注６）。印刷して傍らに置いておくと便利です。
