FrontPage

2011/06/17からのアクセス回数 2934

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/6/

また、Sageのサーバを公開しているサイト(http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/)にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

ベクトルと行列

ベクトルと行列の計算は、幅広い分野で使われています。 Sageを使ってベクトルと行列の特徴的な計算方法について見ていきましょう。

ベクトルの生成

sageでベクトルを生成する関数は、vector関数です。vector関数の使い方を以下に示します。

vector(値のリスト)
または、
vector(環, 値のリスト)

ベクトルの生成例として、$v=(2, 1, 3), w=(1, 1,-4)$の例を以下に示します。環については後ほど紹介します。

sageへの入力:

# ベクトルと行列
v = vector([2,1,3]); v

sageの出力:

(2, 1, 3)

sageへの入力:

w = vector([1,1,-4]); w

sageの出力:

(1, 1, -4)

sageへの入力:

# ベクトルの和
v+w

sageの出力:

(3, 2, -1)

ベクトルの計算

Sageではベクトルの和と差は通常の変数と同じように+, -記号で表すこととができます。

ベクトルのスカラー倍は、ベクトルの各要素にスカラを掛けた形になります。

ベクトルの絶対値は、ベクトル用各要素の自乗の和を平方根となり、ベクトルの距離は2つのベクトルの差の絶対値となります。

このようにSageではベクトルに対して、変数と同じように計算ができる点が特徴です。

sageへの入力:

var('a1 a2 a3 b1 b2 b3 c')
A = vector([a1, a2, a3])
B = vector([b1, b2, b3])
view (A+B)        # ベクトルの和
view (c*A)         # ベクトルのスカラー倍
view (abs(A))      # ベクトルの絶対値
view (abs(A-B))   # ベクトルの距離

out1.png

ベクトルの内積

ベクトルには、内積と外積という特徴的な2つの演算があります。

内積は、2つのベクトルの各要素の積の和として、以下のように定義されます。 $$ \mathbf{v}\cdot\mathbf{w} = \Sigma_{i=1}^N v_i w_i $$

先に生成したベクトルA, Bに対して内積を取ると、定義どおりの結果が得られます。

sageへの入力:

AdB = A.dot_product(B); view(AdB)

out2.png

内積の大きな特徴に、2つのベクトルのなす角度\(\theta\)との間に以下の関係があります。 $$ \mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}| cos \theta $$ この性質を使って、ベクトルの類似を計算するのに、\(cos(\theta)\)がよく使われます。

また、2つのベクトルが直行する場合には、\(cos \theta = 0\)から、 ベクトルの内積はベクトルの直交判定にもよく利用されます。

以下に、\(v=(2, 1, 3), w=(1, 1,-4)\)の内積の結果とベクトルv, wのなす角度\(\theta\)を 求めています。

sageへの入力:

v.dot_product(w)

sageの出力:

-9

sageへの入力:

deg = lambda x : N(x * 180/pi)
# vとwからcos(th)を計算
cos_th = v.dot_product(w)/(abs(v)*abs(w))
th = arccos(cos_th)
print deg(th)

sageの出力:

124.537583786181

3次元プロットで、ベクトルv,wをプロットしたのが以下の図です。 2つのベクトルのなす角度が約120度であることが、見て取れます。

sageへの入力:

ZERO = vector([0, 0, 0])
v_line = line([ZERO, v], rgbcolor='blue')
w_line = line([ZERO, w], rgbcolor='green')
vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red')
(v_line+w_line).show(aspect_ratio=1)

out3.png

ベクトルの外積

同様にベクトルの外積\(\mathbf{v}\times\mathbf{w}\)は、cross_product関数で計算します。

3次元ベクトルAとBの外積の結果は、以下のようなります。

sageへの入力:

AcB = A.cross_product(B); view(AcB)

out4.png

ちょっと覚えにくいので、以下の性質を使って覚えると良いでしょう。 $$ \mathbf{A}\times\mathbf{B} = \left| \begin{array}{ccc} i & j & k \\ a_{1} & a_{2} & a_{3} \\ b_{1} & b_{2} & b_{3} \end{array} \right| $$

x, y, z方向の単位ベクトルU=(i, j, k)を定義し、U, A, Bからなる行列mを作り、 そのdetを求める、結果をそれぞれ、\(i, j, k\)でまとめると、 $$ (a_2 b_3 - a_3 b_2) i + (-a_1 b_3 + a_3 b_1) j + (a_1 b_2 - a_2 b_1) k $$ と先の外積の結果と対応することが分かります。

sageへの入力:

var('i j k')
U = vector([i, j, k])
m = matrix([U, A, B]); view(m)
dm = det(m); view(expand(dm))

out5.png

ベクトルの外積の図形的特徴は、外積の方向はベクトルvからベクトルwにねじを回して進む方向となり、 その大きさはベクトルvとベクトルwの作る平行四辺形の面積となります。 $$ |\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|sin(\theta) $$

Sageを使ってベクトルv, wの外積を求め、その値が \(|v||w| cos (\theta)\)と一致することを見てみましょう。

sageへの入力:

VcW = v.cross_product(w)
# thは内積で求めた結果を利用
print VcW, N(abs(VcW)), N(abs(v)*abs(w)*sin(th))

sageの出力:

(-7, 11, 1) 13.0766968306220 13.0766968306220

3次元プロットで、ベクトルv,wと外積の結果をプロットしたのが以下の図です。 図を少し回転すると、ベクトルvからベクトルwに回したときのねじの進行方向に外積ベクトルがのびている ことが確認できます。

sageへの入力:

vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red')
# 図を回転しないと確認しずらいが、vからwにねじを回した方向になっている
(v_line+w_line+vw_line).show(aspect_ratio=1)

out6.png

行列の生成

行列は、matrix関数で生成します。

matrix(行列の要素のリスト)
または
matrix(環, 行列の要素のリスト)

生成された行列の要素は、配列と同じようにカギ括弧[]に要素のインデックスを指定することで、 参照できます。行を取得する場合には、行のインデックスを指定し、列を取得するにはcolumnメソッド 列のインデックスを指定することで所望の情報を得ることができます。

sageへの入力:

M = matrix([[1,2,3],[3,2,1],[1,2,1]]); view(M)

out7.png

sageへの入力:

print M[1]       # 2行目(インデックスでは1)を取得
print M[0, 2]    # 1行目3列の要素を取得
print M.column(1)# 2列目を取得

sageの出力:

(3, 2, 1)
3
(2, 2, 2)

行列とベクトルの積は、*演算子を使って行い、その結果としてベクトルが返されます。

sageへの入力:

Mw = M*w; print M*w
type(Mw)

sageの出力:

(-9, 1, -1)
<type 'sage.modules.vector_integer_dense.Vector_integer_dense'>

行列の基本演算

行列の基本計算をSageを使ってみてみましょう。行列AとBの和、スカラー倍、積の結果を 以下に示します。行列の積は、順序によって結果が異なることに注意してください。

sageへの入力:

var('a11 a12 a21 a22 b11 b12 b21 b22')
A = matrix([[a11, a12], [a21, a22]])
B = matrix([[b11, b12], [b21, b22]])
view(A)
view(B)
view(A+B)	# 行列の和
view(c*A)	# 行列のスカラー倍
view(A*B)	# 行列の積(AB順)
view(B*A)	# 行列の積(BA順)

out8.png

単位行列

単位行列は、identity_matrixで生成します。identity_matrixの引数には、行列の次数を指定します。

sageへの入力:

Im = identity_matrix(3); view(Im)

out9.png

対角行列

対角行列は、diagonal_matrixで生成します。

diagonal_matrix(対角要素のリスト)

sageへの入力:

D = diagonal_matrix([1, 2, 3]); view(D)

out10.png

行列の解

行列MにベクトルXを掛けて、ベクトルYとなる場合、行列MとベクトルYからXを求めるメソッドがsolve_rightです。 $$ \mathbf{M}\mathbf{X} = \mathbf{Y} $$

solve_rightの例を以下に示します。

sageへの入力:

M = matrix([[1,2,3],[3,2,1],[1,2,1]]); show(M)
Y = vector([0,-4,-1]); show(Y)

out11.png

sageへの入力:

# solve_rightを使って求める
X = M.solve_right(Y); view(X)
# octaveの左除算オペレータ\を使って求める
X = M \ Y; view(X)

out12.png

\(\mathbf{M}\mathbf{X}\)を計算すると、ベクトルY(0, -4, -1)となることから解が正しいことが確認できます。

sageへの入力:

M*X

sageの出力:

(0, -4, -1)

連立方程式と行列の解の関係

先の行列式は、以下のような連立方程式と一致します。 $$ \left\{\begin{array}{rrr} x_1 + 2 x_2 + 3 x_3 & = & 0 \\ 3 x_1 + 2 x_2 + x_3 & = & -4 \\ x_1 + 2 x_2 + x_3 & = & -1 \end{array}\right. $$

上記の連立方程式をsolve関数で計算結果とXの値が同じになることをSageで確かめてみましょう。

sageへの入力:

(x1, x2, x3) = var('x1,x2, x3')
eq = [ [ x1 + 2*x2 + 3*x3 == 0], [3*x1 + 2*x2 + x3 == -4], [x1 + 2*x2 + x3 == -1]]
sol = solve(eq, [x1, x2, x3]); show(sol)

out13.png

転置行列

転置行列は、transpose関数で取得できます。

転置行列の性質は、

  • \((\mathbf{A}^T)^T = \mathbf{A}\)
  • \((\mathbf{A}+\mathbf{B})^T = \mathbf{A}^T + \mathbf{B}^T\)
  • \((\mathbf{A}\mathbf{B})^T = \mathbf{B}^T\mathbf{A}^T\)

最後の性質は、行列の掛け合わせる順序を入れ替えるときに便利です。

sageへの入力:

At = A.transpose(); show(At)           # 転置行列
Bt = B.transpose(); show(Bt)

out14.png

sageへの入力:

view("$(A^T)^T = A$")
show(At.transpose())
view("$(A+B)^T = A^T + B^T$")
show(At + Bt)
show((A+B).transpose())

out15.png

sageへの入力:

view("$(AB)^T = B^T A^T$")
show((A*B).transpose())
show(At*Bt)

out16.png

行列式

行列式detは、以下のように定義されます。 $$ det(\mathbf{A}) = \left| \begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{21} \end{array} \right| $$

行列式は、逆行列の計算に使われるので、よく逆行列が存在するか否かの判別に使用されます。

Sageでは、行列式はdet関数を使って計算されます。

sageへの入力:

view(A.det())

out17.png

sageへの入力:

Mdet = M.det(); print Mdet

sageの出力:

8

逆行列

逆行列は、inverse関数で取得できます。

逆行列の性質は、以下の通りです。

  • \((\mathbf{A}^{-1})^{-1} = \mathbf{A}\)
  • \((\mathbf{A}^{T})^{-1} = (\mathbf{A}^{-1})^T\)
  • \((\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}\)

sageへの入力:

Minv = M.inverse(); view(Minv)

out18.png

sageへの入力:

view("$(A^{-1})^{-1} = A$")
view(Minv.inverse())

out19.png

sageへの入力:

view("$(A^{T})^{-1} = (A^{-1})^{T}$")
view(M.transpose().inverse())
view(M.inverse().transpose())

out20.png

環は、数の集合を表します。 よく使われる環を以下に示します。

  • 整数: ZZ
  • 実数: R
  • 有理数: QQ
  • 複素数: CC
  • 倍精度小数(real double): RDF

例として、有理数の環を使って行列を生成してみます。

sageへの入力:

A = matrix(QQ,3,3,[[2, 4, 0],[3, 1, 0], [0, 1, 1]]); view(A)

out21.png

固有値解析

与えられた正方行列Aの固有値と固有ベクトルは、つぎのように求めます。

行列Aの固有方程式が0ベクトル以外の解を持つときに、 $\lambda$を$\mathbf{x}$の固有値、$x$を固有ベクトルといいます。 $$ \mathbf{A}\mathbf{x} = \lambda\mathbf{x} $$

固有値は、以下の固有方程式の解です。ここで、Eは単位行列です。 $$ det(\mathbf{A} - \lambda\mathbf{E}) = 0 $$

各々の固有値$\lambda_i$に対して、以下の式を満たす固有ベトル$x_i$を求めます。 $$ (\mathbf{A} - \lambda_i \mathbf{E})x_i = 0 $$

例として、以下のような行列Aを使って固有値と固有ベクトルを求めてみます。 $$ \mathbf{A} = \left(\begin{array}{rr} 1 & 3 \\ 2 & 0 \end{array}\right) $$

最初に固有方程式(Mdet)とその解(sol)を求めます。 ここでは、\(\lambda\)の代わりに変数rとして、固有値を求めています。 rの解として、-2と3が求まりました。

sageへの入力:

r = var('r')
E = identity_matrix(2)
A = matrix([[1, 3], [2, 0]])
eq = A - r*E
Mdet = det(eq); view (Mdet.expand())
sol = solve(Mdet, r); view (sol)

out22.png

次に、各固有値に対する固有ベクトルを求めます。方程式eqにr=-2を代入し、 $$ \mathbf{A}_1 \mathbf{x}_1 = 0 $$ を満たす\(\mathbf{x}_1\)を求めるとエラーとなります。 行列A1の1行目と2行目が比例関係であり、2x2の行列A1のランクが1となっているからです。

sageへの入力:

# r = -2に対する固有ベクトル
A1 = eq.subs(r = -2); show(A1)
# A1のrankが1であるため、そのままでは A1*X = Zが求まらない
print A1.rank()
Z = vector([0, 0])
X = A1 \ Z; X

out23.png

そこで、right_kernelメソッドを使って、固有ベクトルを求めると(1, -1)が求まります。

固有ベクトルの定数倍もまた固有ベクトルとなりますので、t(1, -1)のように定数tを付けて 表します。

sageへの入力:

# \の代わりにright_kernelメソッドを使って計算
V = A1.right_kernel(); 
V1 = V.basis_matrix().transpose(); show(V1)
# A1*V1 = 0であることを確認
view(A1*V1)

out24.png

同様に、固有値3に対して、固有ベクトルを求めると、t(1, 2/3)と求まります。

sageへの入力:

# r = 3に対する固有ベクトル
A2 = eq.subs(r = 3); show(A2)
V = A2.right_kernel(); 
V2 = V.basis_matrix().transpose(); show(V2)
# A2*V2 = 0であることを確認
show(A2*V2)

out25.png

行列Aの転置行列を使う理由

sageで固有値と固有ベクトルを取得するには、eigenmatrix_left関数を使います。 eigenmatrix_leftは、与えられた行列Aに対して、以下の関係を満たす対角行列Dと 固有ベクトルP(各固有値に対して、固有ベクトルは行単位で返されることに注意)のタプルを返します。 $$ \mathbf{P} \mathbf{A} = \mathbf{D} \mathbf{P} $$

行列Aの各固有値と固有ベクトルの関係は、以下のようになります。 $$ \begin{array}{ccc} \mathbf{A} \mathbf{x}_1 & = & \lambda_1 \mathbf{x}_1 \\ \mathbf{A} \mathbf{x}_2 & = & \lambda_2 \mathbf{x}_2 \end{array} $$ これを1つにまとめると、以下のようになります。 $$ ( \mathbf{A} \mathbf{x}_1, \mathbf{A} \mathbf{x}_2 ) = ( \lambda_1 \mathbf{x}_1, \lambda_2 \mathbf{x}_2 ) $$ これを行列で表すと、以下のようになります。 $$ \mathbf{A} (\mathbf{x}_1, \mathbf{x}_2) = (\mathbf{x}_1, \mathbf{x}_2) \left(\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_1 \end{array}\right) $$ 両辺の転置行列を取り、\((AB)^T = B^T A^T\)の関係式より、eigenmatrix_left関数と同じ以下の形式を得ます。 $$ \left(\begin{array}{r} \mathbf{x}_1 \\ \mathbf{x}_2 \end{array}\right) \mathbf{A}^T = \left(\begin{array}{cc} \lambda_1 & 0 \\ 0 & \lambda_2 \end{array}\right) \left(\begin{array}{r} \mathbf{x}_1 \\ \mathbf{x}_2 \end{array}\right) $$ 従って、固有ベクトルが固有ベクトルの行単位で返されるeigenmatrix_left関数を使う場合には、 対象となる行列Aの転置行列を使って計算する必要があるのです。

以下に同様の計算をeigenmatrix_left関数を使って求める方法を示します。

sageへの入力:

# 同様の処理をeigenspaces_leftを使って求める
At = A.transpose()
At. eigenspaces_left()

sageの出力:

[
(3, Vector space of degree 2 and dimension 1 over Rational Field
User basis matrix:
[  1 2/3]),
(-2, Vector space of degree 2 and dimension 1 over Rational Field
User basis matrix:
[ 1 -1])
]

計算が正しく求まっていることは、\(\mathbf{P}\mathbf{D}\mathbf{P} = \mathbf{A}^T\) となっていることで確かめることができます。

sageへの入力:

# P*Mt = D*PとなるD, Pを求める
(D, P) = At.eigenmatrix_left()
show(D)
show(P)
# ~P*D*P = Atとなることを確認(~Pは、P.inverse()の省略形)
show(~P*D*P)

out26.png

主成分分析

主成分分析では、固有値の絶対値の大きなものから順に固有ベクトルを求めます。 その結果、固有値の絶対値が小さなものは、行列対する寄与が小さいことから どの固有値までがもとになるAを表現するかを、 固有値に係数を掛けて\(\mathbf{A} = \mathbf{P}^{-1}\mathbf{D}'\mathbf{P}\) 見てみることで主要な成分を求めることができます。

ちょっと乱暴ですが、先の例題の-2の固有値に0.5を掛けて、\(\mathbf{D}'\)を求め、 その影響をみてみましょう。 Aを近似してはいませんが、Aの傾向はつかめていることが分かります。

sageへの入力:

D1 = diagonal_matrix([3, -2*0.5])
show(~P*D1*P)

out27.png

SVD分解

固有値を使った主成分分析は正方行列しか使えないため、n, mが異なる行列での主成分分析にはSVD分解を使います。 $$ \mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^{T} $$ となるU, S, Vを計算するのが、SVD関数です。

以下にSVD関数の例を示します。

sageへの入力:

m = matrix(RDF, 2, range(6)); m
[0.0 1.0 2.0]
[3.0 4.0 5.0]

sageへの入力:

(U, S, V) = m.SVD()
print 'U'
show( U)
print 'S'
show( S)
print 'V'
show( V)

out28.png

計算された\(\mathbf{U}\mathbf{S}\mathbf{V}^{T}\)もかなり精度よく求まっています。

sageへの入力:

show(U*S*V.transpose())

out29.png

コメント

選択肢 投票
おもしろかった 3  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)


添付ファイル: fileout29.png 378件 [詳細] fileout28.png 381件 [詳細] fileout27.png 348件 [詳細] fileout26.png 372件 [詳細] fileout25.png 349件 [詳細] fileout24.png 347件 [詳細] fileout23.png 382件 [詳細] fileout22.png 301件 [詳細] fileout21.png 315件 [詳細] fileout20.png 367件 [詳細] fileout19.png 295件 [詳細] fileout18.png 360件 [詳細] fileout17.png 350件 [詳細] fileout16.png 384件 [詳細] fileout15.png 377件 [詳細] fileout14.png 359件 [詳細] fileout13.png 373件 [詳細] fileout12.png 336件 [詳細] fileout11.png 324件 [詳細] fileout10.png 347件 [詳細] fileout9.png 357件 [詳細] fileout8.png 391件 [詳細] fileout7.png 316件 [詳細] fileout6.png 419件 [詳細] fileout5.png 369件 [詳細] fileout4.png 429件 [詳細] fileout3.png 378件 [詳細] fileout2.png 406件 [詳細] fileout1.png 368件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-04-15 (火) 08:27:08 (1111d)
SmartDoc