CoCalc -- 1049-行列計算に挑戦.sagews

All published worksheets from http://sagenb.org

Project: sagenb.org published worksheets

Path: pub / 1001-1101 / 1049-行列計算に挑戦.sagews

Views: ¹⁶⁸⁷⁰⁴
Image: ubuntu2004

行列計算に挑戦

要素: リスト[要素のインデックス]
指定範囲のサブリスト: リスト[開始インデックス:終了インデックス] $\cdots$ 終了インデックスを含めません。
指定以降のサブリスト: リスト[開始インデックス:] 
先頭から指定までのサブリスト: リスト[:終了インデックス] $\cdots$ 終了インデックスを含めません。

print l1[2]
print l1[1:3]
print l1[2:]
print l1[:2]

3
[2, 3]
[3, 4]
[1, 2]

リストの連結は、+で行います。

l2 = [9, 8, 7]
l1 + l2

[1, 2, 3, 4, 9, 8, 7]

範囲生成関数rangeを使うと簡単に指定した範囲のリストを生成することができます。

range(終了) : ０から終了の前までの整数のリストを生成します
range(開始, 終了) : 開始インデックスから終了インデックスの前まで整数のリストを生成します
range(開始, 終了, 増分) : 開始から終了の前までの値を増分ごとに追加したリストを生成します

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

print range(3)
print range(1, 4)
print range(-1, 7, 2)

[0, 1, 2]
[1, 2, 3]
[-1, 1, 3, 5]

リスト内包は、式とそれにつづくfor節からなり、各要素に対し、式を評価した結果をリストとして返します。

以下の例は、リストl1の各要素の自乗を要素とするリストを返します。

[s^2 for s in l1]

[1, 4, 9, 16]

同様の処理をfor文を使って処理スト以下のようになります。リスト内包が便利なことが分かります。

l3=[]
for s in l1:
    l3.append(s^2)
print l3

[1, 4, 9, 16]

もう一つ、各要素のルートを取る例を示します。

sqrtを式のままリストに保持
sqrtの値を5桁の数値で保持
リストを整形して出力した場合

s1 = [sqrt(s) for s in l1]
print s1
print [n(s, digits=5) for s in s1]
view(s1)

[1, sqrt(2), sqrt(3), 2]
[1.0000, 1.4142, 1.7320, 2.0000]

\newcommand{\Bold}[1]{\mathbf{#1}}\left[1, \sqrt{2}, \sqrt{3}, 2\right]

s1 = [sqrt(s) for s in l1]
print s1
print [n(s, digits=5) for s in s1]
view(s1)

[1, sqrt(2), sqrt(3), 2]
[1.0000, 1.4142, 1.7320, 2.0000]

\newcommand{\Bold}[1]{\mathbf{#1}}\left[1, \sqrt{2}, \sqrt{3}, 2\right]

2つのリストの各要素をリストとして結合したい場合に、zipを使います。

以下の例では、数値と文字のリストをzipで結合します。

t1 = ['a', 'b', 'c', 'd']
print t1
zip(l1, t1)

['a', 'b', 'c', 'd']
[(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd')]

sum関数は、リストの値の和を返します。

# l1の1文字目はLです！
sum(l1)

10

flattenは、入れ子になったリストのすべての要素を1個のリストに並べます。

print [l1, l2]
flatten([l1+l2])

[[1, 2, 3, 4], [9, 8, 7]]
[1, 2, 3, 4, 9, 8, 7]

list_plot([i^2 for i in range(25)], xmax=25)

リストプロット

これで、ようやくリストプロットの例を説明することができます。

list_plot(プロットするリスト)

０から２４までの値の自乗の値をプロットします。

list_plot([i^2 for i in range(25)], xmax=25, plotjoined=true)

ベクトルと行列

ベクトルの生成

ベクトルは、vector関数で生成します。

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

以下は、

v=(2, 1, 3), v=(1, 1,-4)

を生成しています。環については後ほど紹介します。

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

(2, 1, 3)

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

(1, 1, -4)

ベクトルの内積

\mathbf{v}\cdot\mathbf{w}

は、dot_product関数で計算します。

よく知られているようにベクトルの内積は、

$\mathbf{v}\cdot\mathbf{w} = |\mathbf{v}||\mathbf{w}|cos(\theta)$

で表されます。

v.dot_product(w)

-9

ベクトルの外積

同様にベクトルの外積

\mathbf{v}\times\mathbf{w}

は、cross_product関数で計算します。ベクトルの外積の大きさは、

$|\mathbf{v}\times\mathbf{w}| = |\mathbf{v}||\mathbf{w}|sin(\theta)$

であり、その方向は、ベクトルvからベクトルwにねじを回して進む方向となります。

v.cross_product(w)

(-7, 11, 1)

行列の生成

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

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

A = matrix([[1,2,3],[3,2,1],[1,2,1]]); A

[1 2 3]
[3 2 1]
[1 2 1]

行列、ベクトルの積は、*を使って表現します。

w*A

(0, -4, 0)

行列の解

solve_right関数を使って行列の解を得ることができます。

$\mathbf{A}\mathbf{X} = \mathbf{Y}$

となるようなXを求めるのが、solve_right関数です。

Y = vector([0,-4,-1])

X = A.solve_right(Y); X

(-3/2, 0, 1/2)

\mathbf{A}\mathbf{X}

を計算して、解が正しいことを確認しましょう。

A*X

(0, -4, -1)

転置行列

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

最後の、転置行列の積の掛ける順序を変えるための公式には、座標変換行列の処理でお世話になりました。

A.transpose()

[1 3 1]
[2 2 2]
[3 1 1]

行列式

行列式、det関数で取得できます。

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

A.det()

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}

です。

A.inverse()

[   0  1/2 -1/2]
[-1/4 -1/4    1]
[ 1/2    0 -1/2]

単位行列

単位行列は、identity_matrixで生成します。

identity_matrix(3)

[1 0 0]
[0 1 0]
[0 0 1]

対角行列

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

diagonal_matrix(対角要素のリスト)

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

[1 0 0]
[0 2 0]
[0 0 3]

環

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

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

です。以下に有理数の環を使って行列を生成する例を示します。

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

[2 4 0]
[3 1 0]
[0 1 1]

固有値解析

固有値と固有ベクトルは、

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

固有値は、

$det(\mathbf{M} - \lambda\mathbf{E}) = 0$ の解です。

sageで固有値と固有ベクトルを取得するには、eigenmatrix_left関数を使います。eigenmatrix_leftの戻り値は固有値の対角行列と各固有値の固有ベクトルからなる行列Pをタプルとして返します。

先に定義した $\mathbf{A}$ の固有値と固有ベクトルをeigenmatrix_rightで計算してみます。

eigenmatrix_rightは、定義とすこし違いP*self = P*DとなるDとPを返します。

Dの対角成分が固有値で、Pの各列が固有値に対応する固有ベクトルがセットされています。固有値解析の場合、固有値の大きな順に固有ベクトルを計算します。

D, P = A.eigenmatrix_right()
print "D="
print D
print "P="
print P

D=
[ 5  0  0]
[ 0  1  0]
[ 0  0 -2]
P=
[   1    0    1]
[ 3/4    0   -1]
[3/16    1  1/3]

eigenmatrix_rightを確認するために、A*P == D*Pを確認してみます。

A*P == P*D # ~Pは、P.inverse()の省略形

True

主成分分析では、固有値の大きなものを使って処理するため、P*self = D*Pとなる解を返すeigenmatrix_leftの方が、便利なこともあります。

D, P = A.eigenmatrix_left()
print "D="
print D
print "P="
print P

D=
[ 5  0  0]
[ 0  1  0]
[ 0  0 -2]
P=
[   1    1    0]
[   1 -1/3   -4]
[   1 -4/3    0]

元となるAは、

\mathbf{A} = \mathbf{P}^{-1}\mathbf{D}\mathbf{P}

によって表されます。

どの固有値までがもとになるAを表現するかを、固有値に係数を掛けて $\mathbf{A} = \mathbf{P}^{-1}\mathbf{D}'\mathbf{P}$ 見てみることで主要な成分を求めることができます。

A == (~P)*D*P # ~Pは、Pの逆行列

True

求まった固有値の、絶対値の一番小さな1に対して0.5を掛けても

\mathbf{P}^{-1}\mathbf{D}'\mathbf{P}

と

\mathbf{A}

との残差がそんなに大きく変わらないことが分かります。

D=matrix([[5,0,0],[0,1*0.5,0],[0,0,-2]]) # 第2項に0.5を掛けた

print '(~P)*D*P'
print (~P)*D*P
print 'A-(~P)*D*P'
n(A-(~P)*D*P, digits=5)

(~P)*D*P
[ 2.00000000000000  4.00000000000000 0.000000000000000]
[ 3.00000000000000  1.00000000000000 0.000000000000000]
[0.125000000000000 0.958333333333333 0.500000000000000]
A-(~P)*D*P
[4.4409e-16 4.4409e-16    0.00000]
[   0.00000    0.00000    0.00000]
[  -0.12500   0.041667    0.50000]

SVD分解

固有値を使った主成分分析は正方行列しか使えないため、n, mが異なる行列での主成分分析にはSVD分解を使います。

$\mathbf{A} = \mathbf{U}\mathbf{S}\mathbf{V}^{T}$

となるU, S, Vを計算するのが、SVD関数です。

m = matrix(RDF, 2, range(6)); m

[0.0 1.0 2.0]
[3.0 4.0 5.0]

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

U
[-0.274721127897 -0.961523947641]
[-0.961523947641  0.274721127897]
S
[7.34846922835           0.0           0.0]
[          0.0           1.0           0.0]
V
[-0.392540507864  0.824163383692  0.408248290464]
[-0.560772154092  0.137360563949 -0.816496580928]
[ -0.72900380032 -0.549442255795  0.408248290464]

計算された

\mathbf{U}\mathbf{S}\mathbf{V}^{T}

もかなり精度よく求まっています。

U*S*V.transpose()

[9.99200722163e-16               1.0               2.0]
[              3.0               4.0               5.0]

S*V.transpose()

[ -2.88457184292  -4.12081691846    -5.357061994]
[ 0.824163383692  0.137360563949 -0.549442255795]

SS = matrix(2, 3, [7.34846922835, 0, 0, 0, 1*0.5, 0]); SS

[ 7.34846922835000 0.000000000000000 0.000000000000000]
[0.000000000000000 0.500000000000000 0.000000000000000]

print "A'"
print n(U*SS*V.transpose(), digits=5)
print '残差'
print n(m-U*SS*V.transpose(), digits=5)

A'
[0.39623  1.0660  1.7358]
[ 2.8868  3.9811  5.0755]
残差
[ -0.39623 -0.066038   0.26415]
[  0.11321  0.018868 -0.075472]