sage/text/ペクトルと行列のバックアップ(No.2)

バックアップ一覧
差分を表示
現在との差分を表示
ソースを表示
sage/text/ペクトルと行列へ行く。
- 1 (2012-06-17 (日) 17:20:57)
- 2 (2012-06-17 (日) 17:20:57)
- 3 (2013-06-13 (木) 14:21:27)
- 4 (2014-04-13 (日) 02:35:44)

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

ここで紹介した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ではベクトルの和と差は通常の変数と同じように+, -記号で表すこととができます。

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

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

このように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))   # ベクトルの距離

↑

ベクトルの内積 †

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

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

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

sageへの入力:

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

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

また、２つのベクトルが直行する場合には、$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

３次元プロットで、ベクトルv,wをプロットしたのが以下の図です。２つのベクトルのなす角度が約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)

↑

ベクトルの外積 †

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

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

sageへの入力:

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

ちょっと覚えにくいので、以下の性質を使って覚えると良いでしょう。 $$ \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))

ベクトルの外積の図形的特徴は、外積の方向はベクトル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

３次元プロットで、ベクトル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)

↑

行列の生成 †

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

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

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

sageへの入力:

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

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順）

↑

単位行列 †

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

sageへの入力:

Im = identity_matrix(3); view(Im)

↑

対角行列 †

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

diagonal_matrix(対角要素のリスト)

sageへの入力:

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

↑

行列の解 †

行列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)

sageへの入力:

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

$\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)

↑

転置行列 †

転置行列は、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)

sageへの入力:

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

sageへの入力:

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

↑

行列式 †

行列式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())

sageへの入力:

Mdet = M.det(); print Mdet

sageの出力:

↑

逆行列 †

逆行列は、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)

sageへの入力:

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

sageへの入力:

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

↑

環 †

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

整数: ZZ

実数: R

有理数: QQ

複素数: CC

倍精度小数(real double): RDF

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

sageへの入力:

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

↑

固有値解析 †

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

行列Aの固有方程式が０ベクトル以外の解を持つときに、 $\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と３が求まりました。

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)

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

sageへの入力:

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

そこで、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)

同様に、固有値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)

↑

行列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)

↑

主成分分析 †

主成分分析では、固有値の絶対値の大きなものから順に固有ベクトルを求めます。その結果、固有値の絶対値が小さなものは、行列対する寄与が小さいことからどの固有値までがもとになる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)

↑

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)

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

sageへの入力:

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

↑

コメント †

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

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

sage/text/ペクトルと行列 のバックアップ(No.2)

ベクトルと行列 †

ベクトルの生成 †

ベクトルの計算 †

ベクトルの内積 †

ベクトルの外積 †

行列の生成 †

行列の基本演算 †

単位行列 †

対角行列 †

行列の解 †

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

転置行列 †

行列式 †

逆行列 †

環 †

固有値解析 †

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

主成分分析 †

SVD分解 †

コメント †

sage/text/ペクトルと行列のバックアップ(No.2)