sage/text/ペクトルと行列
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2011/06/17からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/6/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
** ベクトルと行列 [#z648c7a2]
ベクトルと行列の計算は、幅広い分野で使われています。
Sageを使ってベクトルと行列の特徴的な計算方法について見て...
*** ベクトルの生成 [#x8f501d1]
sageでベクトルを生成する関数は、vector関数です。vector関...
#pre{{
vector(値のリスト)
または、
vector(環, 値のリスト)
}}
ベクトルの生成例として、$v=(2, 1, 3), w=(1, 1,-4)$の例を...
sageへの入力:
#pre{{
# ベクトルと行列
v = vector([2,1,3]); v
}}
sageの出力:
#pre{{
(2, 1, 3)
}}
sageへの入力:
#pre{{
w = vector([1,1,-4]); w
}}
sageの出力:
#pre{{
(1, 1, -4)
}}
sageへの入力:
#pre{{
# ベクトルの和
v+w
}}
sageの出力:
#pre{{
(3, 2, -1)
}}
*** ベクトルの計算 [#i27629fe]
Sageではベクトルの和と差は通常の変数と同じように+, -記号...
ベクトルのスカラー倍は、ベクトルの各要素にスカラを掛けた...
ベクトルの絶対値は、ベクトル用各要素の自乗の和を平方根と...
このようにSageではベクトルに対して、変数と同じように計算...
sageへの入力:
#pre{{
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)) # ベクトルの距離
}}
&ref(out1.png);
*** ベクトルの内積 [#l6c177e9]
ベクトルには、内積と外積という特徴的な２つの演算がありま...
内積は、２つのベクトルの各要素の積の和として、以下のよう...
$$
\mathbf{v}\cdot\mathbf{w} = \Sigma_{i=1}^N v_i w_i
$$
先に生成したベクトルA, Bに対して内積を取ると、定義どおり...
sageへの入力:
#pre{{
AdB = A.dot_product(B); view(AdB)
}}
&ref(out2.png);
内積の大きな特徴に、２つのベクトルのなす角度\(\theta\)と...
$$
\mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}| cos ...
$$
この性質を使って、ベクトルの類似を計算するのに、\(cos(\th...
また、２つのベクトルが直行する場合には、\(cos \theta = 0\...
ベクトルの内積はベクトルの直交判定にもよく利用されます。
以下に、\(v=(2, 1, 3), w=(1, 1,-4)\)の内積の結果とベクト...
求めています。
sageへの入力:
#pre{{
v.dot_product(w)
}}
sageの出力:
#pre{{
-9
}}
sageへの入力:
#pre{{
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の出力:
#pre{{
124.537583786181
}}
３次元プロットで、ベクトルv,wをプロットしたのが以下の図で...
２つのベクトルのなす角度が約120度であることが、見て取れま...
sageへの入力:
#pre{{
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)
}}
&ref(out3.png);
*** ベクトルの外積 [#v5e88e45]
同様にベクトルの外積\(\mathbf{v}\times\mathbf{w}\)は、cro...
３次元ベクトルAとBの外積の結果は、以下のようなります。
sageへの入力:
#pre{{
AcB = A.cross_product(B); view(AcB)
}}
&ref(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から...
その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...
$$
と先の外積の結果と対応することが分かります。
sageへの入力:
#pre{{
var('i j k')
U = vector([i, j, k])
m = matrix([U, A, B]); view(m)
dm = det(m); view(expand(dm))
}}
&ref(out5.png);
ベクトルの外積の図形的特徴は、外積の方向はベクトルvからベ...
その大きさはベクトルvとベクトルwの作る平行四辺形の面積と...
$$
|\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|si...
$$
Sageを使ってベクトルv, wの外積を求め、その値が \(|v||w| c...
sageへの入力:
#pre{{
VcW = v.cross_product(w)
# thは内積で求めた結果を利用
print VcW, N(abs(VcW)), N(abs(v)*abs(w)*sin(th))
}}
sageの出力:
#pre{{
(-7, 11, 1) 13.0766968306220 13.0766968306220
}}
３次元プロットで、ベクトルv,wと外積の結果をプロットしたの...
図を少し回転すると、ベクトルvからベクトルwに回したときの...
ことが確認できます。
sageへの入力:
#pre{{
vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red')
# 図を回転しないと確認しずらいが、vからwにねじを回した方...
(v_line+w_line+vw_line).show(aspect_ratio=1)
}}
&ref(out6.png);
*** 行列の生成 [#s3e97f73]
行列は、matrix関数で生成します。
#pre{{
matrix(行列の要素のリスト)
または
matrix(環, 行列の要素のリスト)
}}
生成された行列の要素は、配列と同じようにカギ括弧[]に要素...
参照できます。行を取得する場合には、行のインデックスを指...
列のインデックスを指定することで所望の情報を得ることがで...
sageへの入力:
#pre{{
M = matrix([[1,2,3],[3,2,1],[1,2,1]]); view(M)
}}
&ref(out7.png);
sageへの入力:
#pre{{
print M[1] # 2行目（インデックスでは1）を取得
print M[0, 2] # 1行目3列の要素を取得
print M.column(1)# 2列目を取得
}}
sageの出力:
#pre{{
(3, 2, 1)
3
(2, 2, 2)
}}
行列とベクトルの積は、*演算子を使って行い、その結果として...
sageへの入力:
#pre{{
Mw = M*w; print M*w
type(Mw)
}}
sageの出力:
#pre{{
(-9, 1, -1)
<type 'sage.modules.vector_integer_dense.Vector_integer_d...
}}
*** 行列の基本演算 [#x9cdf52f]
行列の基本計算をSageを使ってみてみましょう。行列AとBの和...
以下に示します。行列の積は、順序によって結果が異なること...
sageへの入力:
#pre{{
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順）
}}
&ref(out8.png);
*** 単位行列 [#s4e1b2d5]
単位行列は、identity_matrixで生成します。identity_matrix...
sageへの入力:
#pre{{
Im = identity_matrix(3); view(Im)
}}
&ref(out9.png);
*** 対角行列 [#hcbe3560]
対角行列は、diagonal_matrixで生成します。
#pre{{
diagonal_matrix(対角要素のリスト)
}}
sageへの入力:
#pre{{
D = diagonal_matrix([1, 2, 3]); view(D)
}}
&ref(out10.png);
*** 行列の解 [#r443124b]
行列MにベクトルXを掛けて、ベクトルYとなる場合、行列Mとベ...
$$
\mathbf{M}\mathbf{X} = \mathbf{Y}
$$
solve_rightの例を以下に示します。
sageへの入力:
#pre{{
M = matrix([[1,2,3],[3,2,1],[1,2,1]]); show(M)
Y = vector([0,-4,-1]); show(Y)
}}
&ref(out11.png);
sageへの入力:
#pre{{
# solve_rightを使って求める
X = M.solve_right(Y); view(X)
# octaveの左除算オペレータ\を使って求める
X = M \ Y; view(X)
}}
&ref(out12.png);
\(\mathbf{M}\mathbf{X}\)を計算すると、ベクトルY(0, -4, -1...
sageへの入力:
#pre{{
M*X
}}
sageの出力:
#pre{{
(0, -4, -1)
}}
*** 連立方程式と行列の解の関係 [#g788f4e4]
先の行列式は、以下のような連立方程式と一致します。
$$
\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への入力:
#pre{{
(x1, x2, x3) = var('x1,x2, x3')
eq = [ [ x1 + 2*x2 + 3*x3 == 0], [3*x1 + 2*x2 + x3 == -4]...
sol = solve(eq, [x1, x2, x3]); show(sol)
}}
&ref(out13.png);
*** 転置行列 [#c15e6cae]
転置行列は、transpose関数で取得できます。
転置行列の性質は、
- \((\mathbf{A}^T)^T = \mathbf{A}\)
- \((\mathbf{A}+\mathbf{B})^T = \mathbf{A}^T + \mathbf{B}...
- \((\mathbf{A}\mathbf{B})^T = \mathbf{B}^T\mathbf{A}^T\)
最後の性質は、行列の掛け合わせる順序を入れ替えるときに便...
sageへの入力:
#pre{{
At = A.transpose(); show(At) # 転置行列
Bt = B.transpose(); show(Bt)
}}
&ref(out14.png);
sageへの入力:
#pre{{
view("$(A^T)^T = A$")
show(At.transpose())
view("$(A+B)^T = A^T + B^T$")
show(At + Bt)
show((A+B).transpose())
}}
&ref(out15.png);
sageへの入力:
#pre{{
view("$(AB)^T = B^T A^T$")
show((A*B).transpose())
show(At*Bt)
}}
&ref(out16.png);
*** 行列式 [#c131b384]
行列式detは、以下のように定義されます。
$$
det(\mathbf{A}) =
\left| \begin{array}{cc}
a_{11} & a_{12} \\
a_{21} & a_{21}
\end{array} \right|
$$
行列式は、逆行列の計算に使われるので、よく逆行列が存在す...
Sageでは、行列式はdet関数を使って計算されます。
sageへの入力:
#pre{{
view(A.det())
}}
&ref(out17.png);
sageへの入力:
#pre{{
Mdet = M.det(); print Mdet
}}
sageの出力:
#pre{{
8
}}
*** 逆行列 [#r64cbdb6]
逆行列は、inverse関数で取得できます。
逆行列の性質は、以下の通りです。
- \((\mathbf{A}^{-1})^{-1} = \mathbf{A}\)
- \((\mathbf{A}^{T})^{-1} = (\mathbf{A}^{-1})^T\)
- \((\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{...
sageへの入力:
#pre{{
Minv = M.inverse(); view(Minv)
}}
&ref(out18.png);
sageへの入力:
#pre{{
view("$(A^{-1})^{-1} = A$")
view(Minv.inverse())
}}
&ref(out19.png);
sageへの入力:
#pre{{
view("$(A^{T})^{-1} = (A^{-1})^{T}$")
view(M.transpose().inverse())
view(M.inverse().transpose())
}}
&ref(out20.png);
*** 環 [#o677bcd5]
環は、数の集合を表します。
よく使われる環を以下に示します。
- 整数: ZZ
- 実数: R
- 有理数: QQ
- 複素数: CC
- 倍精度小数(real double): RDF
例として、有理数の環を使って行列を生成してみます。
sageへの入力:
#pre{{
A = matrix(QQ,3,3,[[2, 4, 0],[3, 1, 0], [0, 1, 1]]); view...
}}
&ref(out21.png);
*** 固有値解析 [#ob68f797]
与えられた正方行列Aの固有値と固有ベクトルは、つぎのように...
行列Aの固有方程式が０ベクトル以外の解を持つときに、
$\lambda$を$\mathbf{x}$の固有値、$x$を固有ベクトルといい...
$$
\mathbf{A}\mathbf{x} = \lambda\mathbf{x}
$$
固有値は、以下の固有方程式の解です。ここで、Eは単位行列で...
$$
det(\mathbf{A} - \lambda\mathbf{E}) = 0
$$
各々の固有値$\lambda_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への入力:
#pre{{
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)
}}
&ref(out22.png);
次に、各固有値に対する固有ベクトルを求めます。方程式eqにr...
$$
\mathbf{A}_1 \mathbf{x}_1 = 0
$$
を満たす\(\mathbf{x}_1\)を求めるとエラーとなります。
行列A1の１行目と２行目が比例関係であり、2x2の行列A1のラン...
sageへの入力:
#pre{{
# 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
}}
&ref(out23.png);
そこで、right_kernelメソッドを使って、固有ベクトルを求め...
固有ベクトルの定数倍もまた固有ベクトルとなりますので、t(1...
表します。
sageへの入力:
#pre{{
# \の代わりにright_kernelメソッドを使って計算
V = A1.right_kernel();
V1 = V.basis_matrix().transpose(); show(V1)
# A1*V1 = 0であることを確認
view(A1*V1)
}}
&ref(out24.png);
同様に、固有値3に対して、固有ベクトルを求めると、t(1, 2/3...
sageへの入力:
#pre{{
# 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)
}}
&ref(out25.png);
*** 行列Aの転置行列を使う理由 [#g7fc3ed9]
sageで固有値と固有ベクトルを取得するには、eigenmatrix_lef...
eigenmatrix_leftは、与えられた行列Aに対して、以下の関係を...
固有ベクトル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 ) = ( ...
$$
これを行列で表すと、以下のようになります。
$$
\mathbf{A} (\mathbf{x}_1, \mathbf{x}_2) = (\mathbf{x}_1, ...
\left(\begin{array}{cc}
\lambda_1 & 0 \\
0 & \lambda_1
\end{array}\right)
$$
両辺の転置行列を取り、\((AB)^T = B^T A^T\)の関係式より、e...
$$
\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)
$$
従って、固有ベクトルが固有ベクトルの行単位で返されるeigen...
対象となる行列Aの転置行列を使って計算する必要があるのです。
以下に同様の計算をeigenmatrix_left関数を使って求める方法...
sageへの入力:
#pre{{
# 同様の処理をeigenspaces_leftを使って求める
At = A.transpose()
At. eigenspaces_left()
}}
sageの出力:
#pre{{
[
(3, Vector space of degree 2 and dimension 1 over Rationa...
User basis matrix:
[ 1 2/3]),
(-2, Vector space of degree 2 and dimension 1 over Ration...
User basis matrix:
[ 1 -1])
]
}}
計算が正しく求まっていることは、\(\mathbf{P}\mathbf{D}\ma...
となっていることで確かめることができます。
sageへの入力:
#pre{{
# 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)
}}
&ref(out26.png);
*** 主成分分析 [#w0cd92de]
主成分分析では、固有値の絶対値の大きなものから順に固有ベ...
その結果、固有値の絶対値が小さなものは、行列対する寄与が...
どの固有値までがもとになるAを表現するかを、
固有値に係数を掛けて\(\mathbf{A} = \mathbf{P}^{-1}\mathbf...
見てみることで主要な成分を求めることができます。
ちょっと乱暴ですが、先の例題の-2の固有値に0.5を掛けて、\(...
その影響をみてみましょう。
Aを近似してはいませんが、Aの傾向はつかめていることが分か...
sageへの入力:
#pre{{
D1 = diagonal_matrix([3, -2*0.5])
show(~P*D1*P)
}}
&ref(out27.png);
*** SVD分解 [#b390fce5]
固有値を使った主成分分析は正方行列しか使えないため、n, m...
$$
\mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^{T}
$$
となるU, S, Vを計算するのが、SVD関数です。
以下にSVD関数の例を示します。
sageへの入力:
#pre{{
m = matrix(RDF, 2, range(6)); m
}}
#pre{{
[0.0 1.0 2.0]
[3.0 4.0 5.0]
}}
sageへの入力:
#pre{{
(U, S, V) = m.SVD()
print 'U'
show( U)
print 'S'
show( S)
print 'V'
show( V)
}}
&ref(out28.png);
計算された\(\mathbf{U}\mathbf{S}\mathbf{V}^{T}\)もかなり...
sageへの入力:
#pre{{
show(U*S*V.transpose())
}}
&ref(out29.png);
** コメント [#e9aaf72c]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2011/06/17からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/6/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
** ベクトルと行列 [#z648c7a2]
ベクトルと行列の計算は、幅広い分野で使われています。
Sageを使ってベクトルと行列の特徴的な計算方法について見て...
*** ベクトルの生成 [#x8f501d1]
sageでベクトルを生成する関数は、vector関数です。vector関...
#pre{{
vector(値のリスト)
または、
vector(環, 値のリスト)
}}
ベクトルの生成例として、$v=(2, 1, 3), w=(1, 1,-4)$の例を...
sageへの入力:
#pre{{
# ベクトルと行列
v = vector([2,1,3]); v
}}
sageの出力:
#pre{{
(2, 1, 3)
}}
sageへの入力:
#pre{{
w = vector([1,1,-4]); w
}}
sageの出力:
#pre{{
(1, 1, -4)
}}
sageへの入力:
#pre{{
# ベクトルの和
v+w
}}
sageの出力:
#pre{{
(3, 2, -1)
}}
*** ベクトルの計算 [#i27629fe]
Sageではベクトルの和と差は通常の変数と同じように+, -記号...
ベクトルのスカラー倍は、ベクトルの各要素にスカラを掛けた...
ベクトルの絶対値は、ベクトル用各要素の自乗の和を平方根と...
このようにSageではベクトルに対して、変数と同じように計算...
sageへの入力:
#pre{{
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)) # ベクトルの距離
}}
&ref(out1.png);
*** ベクトルの内積 [#l6c177e9]
ベクトルには、内積と外積という特徴的な２つの演算がありま...
内積は、２つのベクトルの各要素の積の和として、以下のよう...
$$
\mathbf{v}\cdot\mathbf{w} = \Sigma_{i=1}^N v_i w_i
$$
先に生成したベクトルA, Bに対して内積を取ると、定義どおり...
sageへの入力:
#pre{{
AdB = A.dot_product(B); view(AdB)
}}
&ref(out2.png);
内積の大きな特徴に、２つのベクトルのなす角度\(\theta\)と...
$$
\mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}| cos ...
$$
この性質を使って、ベクトルの類似を計算するのに、\(cos(\th...
また、２つのベクトルが直行する場合には、\(cos \theta = 0\...
ベクトルの内積はベクトルの直交判定にもよく利用されます。
以下に、\(v=(2, 1, 3), w=(1, 1,-4)\)の内積の結果とベクト...
求めています。
sageへの入力:
#pre{{
v.dot_product(w)
}}
sageの出力:
#pre{{
-9
}}
sageへの入力:
#pre{{
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の出力:
#pre{{
124.537583786181
}}
３次元プロットで、ベクトルv,wをプロットしたのが以下の図で...
２つのベクトルのなす角度が約120度であることが、見て取れま...
sageへの入力:
#pre{{
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)
}}
&ref(out3.png);
*** ベクトルの外積 [#v5e88e45]
同様にベクトルの外積\(\mathbf{v}\times\mathbf{w}\)は、cro...
３次元ベクトルAとBの外積の結果は、以下のようなります。
sageへの入力:
#pre{{
AcB = A.cross_product(B); view(AcB)
}}
&ref(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から...
その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...
$$
と先の外積の結果と対応することが分かります。
sageへの入力:
#pre{{
var('i j k')
U = vector([i, j, k])
m = matrix([U, A, B]); view(m)
dm = det(m); view(expand(dm))
}}
&ref(out5.png);
ベクトルの外積の図形的特徴は、外積の方向はベクトルvからベ...
その大きさはベクトルvとベクトルwの作る平行四辺形の面積と...
$$
|\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|si...
$$
Sageを使ってベクトルv, wの外積を求め、その値が \(|v||w| c...
sageへの入力:
#pre{{
VcW = v.cross_product(w)
# thは内積で求めた結果を利用
print VcW, N(abs(VcW)), N(abs(v)*abs(w)*sin(th))
}}
sageの出力:
#pre{{
(-7, 11, 1) 13.0766968306220 13.0766968306220
}}
３次元プロットで、ベクトルv,wと外積の結果をプロットしたの...
図を少し回転すると、ベクトルvからベクトルwに回したときの...
ことが確認できます。
sageへの入力:
#pre{{
vw_line = line([ZERO, v.cross_product(w)], rgbcolor='red')
# 図を回転しないと確認しずらいが、vからwにねじを回した方...
(v_line+w_line+vw_line).show(aspect_ratio=1)
}}
&ref(out6.png);
*** 行列の生成 [#s3e97f73]
行列は、matrix関数で生成します。
#pre{{
matrix(行列の要素のリスト)
または
matrix(環, 行列の要素のリスト)
}}
生成された行列の要素は、配列と同じようにカギ括弧[]に要素...
参照できます。行を取得する場合には、行のインデックスを指...
列のインデックスを指定することで所望の情報を得ることがで...
sageへの入力:
#pre{{
M = matrix([[1,2,3],[3,2,1],[1,2,1]]); view(M)
}}
&ref(out7.png);
sageへの入力:
#pre{{
print M[1] # 2行目（インデックスでは1）を取得
print M[0, 2] # 1行目3列の要素を取得
print M.column(1)# 2列目を取得
}}
sageの出力:
#pre{{
(3, 2, 1)
3
(2, 2, 2)
}}
行列とベクトルの積は、*演算子を使って行い、その結果として...
sageへの入力:
#pre{{
Mw = M*w; print M*w
type(Mw)
}}
sageの出力:
#pre{{
(-9, 1, -1)
<type 'sage.modules.vector_integer_dense.Vector_integer_d...
}}
*** 行列の基本演算 [#x9cdf52f]
行列の基本計算をSageを使ってみてみましょう。行列AとBの和...
以下に示します。行列の積は、順序によって結果が異なること...
sageへの入力:
#pre{{
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順）
}}
&ref(out8.png);
*** 単位行列 [#s4e1b2d5]
単位行列は、identity_matrixで生成します。identity_matrix...
sageへの入力:
#pre{{
Im = identity_matrix(3); view(Im)
}}
&ref(out9.png);
*** 対角行列 [#hcbe3560]
対角行列は、diagonal_matrixで生成します。
#pre{{
diagonal_matrix(対角要素のリスト)
}}
sageへの入力:
#pre{{
D = diagonal_matrix([1, 2, 3]); view(D)
}}
&ref(out10.png);
*** 行列の解 [#r443124b]
行列MにベクトルXを掛けて、ベクトルYとなる場合、行列Mとベ...
$$
\mathbf{M}\mathbf{X} = \mathbf{Y}
$$
solve_rightの例を以下に示します。
sageへの入力:
#pre{{
M = matrix([[1,2,3],[3,2,1],[1,2,1]]); show(M)
Y = vector([0,-4,-1]); show(Y)
}}
&ref(out11.png);
sageへの入力:
#pre{{
# solve_rightを使って求める
X = M.solve_right(Y); view(X)
# octaveの左除算オペレータ\を使って求める
X = M \ Y; view(X)
}}
&ref(out12.png);
\(\mathbf{M}\mathbf{X}\)を計算すると、ベクトルY(0, -4, -1...
sageへの入力:
#pre{{
M*X
}}
sageの出力:
#pre{{
(0, -4, -1)
}}
*** 連立方程式と行列の解の関係 [#g788f4e4]
先の行列式は、以下のような連立方程式と一致します。
$$
\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への入力:
#pre{{
(x1, x2, x3) = var('x1,x2, x3')
eq = [ [ x1 + 2*x2 + 3*x3 == 0], [3*x1 + 2*x2 + x3 == -4]...
sol = solve(eq, [x1, x2, x3]); show(sol)
}}
&ref(out13.png);
*** 転置行列 [#c15e6cae]
転置行列は、transpose関数で取得できます。
転置行列の性質は、
- \((\mathbf{A}^T)^T = \mathbf{A}\)
- \((\mathbf{A}+\mathbf{B})^T = \mathbf{A}^T + \mathbf{B}...
- \((\mathbf{A}\mathbf{B})^T = \mathbf{B}^T\mathbf{A}^T\)
最後の性質は、行列の掛け合わせる順序を入れ替えるときに便...
sageへの入力:
#pre{{
At = A.transpose(); show(At) # 転置行列
Bt = B.transpose(); show(Bt)
}}
&ref(out14.png);
sageへの入力:
#pre{{
view("$(A^T)^T = A$")
show(At.transpose())
view("$(A+B)^T = A^T + B^T$")
show(At + Bt)
show((A+B).transpose())
}}
&ref(out15.png);
sageへの入力:
#pre{{
view("$(AB)^T = B^T A^T$")
show((A*B).transpose())
show(At*Bt)
}}
&ref(out16.png);
*** 行列式 [#c131b384]
行列式detは、以下のように定義されます。
$$
det(\mathbf{A}) =
\left| \begin{array}{cc}
a_{11} & a_{12} \\
a_{21} & a_{21}
\end{array} \right|
$$
行列式は、逆行列の計算に使われるので、よく逆行列が存在す...
Sageでは、行列式はdet関数を使って計算されます。
sageへの入力:
#pre{{
view(A.det())
}}
&ref(out17.png);
sageへの入力:
#pre{{
Mdet = M.det(); print Mdet
}}
sageの出力:
#pre{{
8
}}
*** 逆行列 [#r64cbdb6]
逆行列は、inverse関数で取得できます。
逆行列の性質は、以下の通りです。
- \((\mathbf{A}^{-1})^{-1} = \mathbf{A}\)
- \((\mathbf{A}^{T})^{-1} = (\mathbf{A}^{-1})^T\)
- \((\mathbf{A}\mathbf{B})^{-1} = \mathbf{B}^{-1}\mathbf{...
sageへの入力:
#pre{{
Minv = M.inverse(); view(Minv)
}}
&ref(out18.png);
sageへの入力:
#pre{{
view("$(A^{-1})^{-1} = A$")
view(Minv.inverse())
}}
&ref(out19.png);
sageへの入力:
#pre{{
view("$(A^{T})^{-1} = (A^{-1})^{T}$")
view(M.transpose().inverse())
view(M.inverse().transpose())
}}
&ref(out20.png);
*** 環 [#o677bcd5]
環は、数の集合を表します。
よく使われる環を以下に示します。
- 整数: ZZ
- 実数: R
- 有理数: QQ
- 複素数: CC
- 倍精度小数(real double): RDF
例として、有理数の環を使って行列を生成してみます。
sageへの入力:
#pre{{
A = matrix(QQ,3,3,[[2, 4, 0],[3, 1, 0], [0, 1, 1]]); view...
}}
&ref(out21.png);
*** 固有値解析 [#ob68f797]
与えられた正方行列Aの固有値と固有ベクトルは、つぎのように...
行列Aの固有方程式が０ベクトル以外の解を持つときに、
$\lambda$を$\mathbf{x}$の固有値、$x$を固有ベクトルといい...
$$
\mathbf{A}\mathbf{x} = \lambda\mathbf{x}
$$
固有値は、以下の固有方程式の解です。ここで、Eは単位行列で...
$$
det(\mathbf{A} - \lambda\mathbf{E}) = 0
$$
各々の固有値$\lambda_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への入力:
#pre{{
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)
}}
&ref(out22.png);
次に、各固有値に対する固有ベクトルを求めます。方程式eqにr...
$$
\mathbf{A}_1 \mathbf{x}_1 = 0
$$
を満たす\(\mathbf{x}_1\)を求めるとエラーとなります。
行列A1の１行目と２行目が比例関係であり、2x2の行列A1のラン...
sageへの入力:
#pre{{
# 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
}}
&ref(out23.png);
そこで、right_kernelメソッドを使って、固有ベクトルを求め...
固有ベクトルの定数倍もまた固有ベクトルとなりますので、t(1...
表します。
sageへの入力:
#pre{{
# \の代わりにright_kernelメソッドを使って計算
V = A1.right_kernel();
V1 = V.basis_matrix().transpose(); show(V1)
# A1*V1 = 0であることを確認
view(A1*V1)
}}
&ref(out24.png);
同様に、固有値3に対して、固有ベクトルを求めると、t(1, 2/3...
sageへの入力:
#pre{{
# 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)
}}
&ref(out25.png);
*** 行列Aの転置行列を使う理由 [#g7fc3ed9]
sageで固有値と固有ベクトルを取得するには、eigenmatrix_lef...
eigenmatrix_leftは、与えられた行列Aに対して、以下の関係を...
固有ベクトル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 ) = ( ...
$$
これを行列で表すと、以下のようになります。
$$
\mathbf{A} (\mathbf{x}_1, \mathbf{x}_2) = (\mathbf{x}_1, ...
\left(\begin{array}{cc}
\lambda_1 & 0 \\
0 & \lambda_1
\end{array}\right)
$$
両辺の転置行列を取り、\((AB)^T = B^T A^T\)の関係式より、e...
$$
\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)
$$
従って、固有ベクトルが固有ベクトルの行単位で返されるeigen...
対象となる行列Aの転置行列を使って計算する必要があるのです。
以下に同様の計算をeigenmatrix_left関数を使って求める方法...
sageへの入力:
#pre{{
# 同様の処理をeigenspaces_leftを使って求める
At = A.transpose()
At. eigenspaces_left()
}}
sageの出力:
#pre{{
[
(3, Vector space of degree 2 and dimension 1 over Rationa...
User basis matrix:
[ 1 2/3]),
(-2, Vector space of degree 2 and dimension 1 over Ration...
User basis matrix:
[ 1 -1])
]
}}
計算が正しく求まっていることは、\(\mathbf{P}\mathbf{D}\ma...
となっていることで確かめることができます。
sageへの入力:
#pre{{
# 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)
}}
&ref(out26.png);
*** 主成分分析 [#w0cd92de]
主成分分析では、固有値の絶対値の大きなものから順に固有ベ...
その結果、固有値の絶対値が小さなものは、行列対する寄与が...
どの固有値までがもとになるAを表現するかを、
固有値に係数を掛けて\(\mathbf{A} = \mathbf{P}^{-1}\mathbf...
見てみることで主要な成分を求めることができます。
ちょっと乱暴ですが、先の例題の-2の固有値に0.5を掛けて、\(...
その影響をみてみましょう。
Aを近似してはいませんが、Aの傾向はつかめていることが分か...
sageへの入力:
#pre{{
D1 = diagonal_matrix([3, -2*0.5])
show(~P*D1*P)
}}
&ref(out27.png);
*** SVD分解 [#b390fce5]
固有値を使った主成分分析は正方行列しか使えないため、n, m...
$$
\mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^{T}
$$
となるU, S, Vを計算するのが、SVD関数です。
以下にSVD関数の例を示します。
sageへの入力:
#pre{{
m = matrix(RDF, 2, range(6)); m
}}
#pre{{
[0.0 1.0 2.0]
[3.0 4.0 5.0]
}}
sageへの入力:
#pre{{
(U, S, V) = m.SVD()
print 'U'
show( U)
print 'S'
show( S)
print 'V'
show( V)
}}
&ref(out28.png);
計算された\(\mathbf{U}\mathbf{S}\mathbf{V}^{T}\)もかなり...
sageへの入力:
#pre{{
show(U*S*V.transpose())
}}
&ref(out29.png);
** コメント [#e9aaf72c]
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
