線形代数(Linear algebra)

cc by Shigeto R. Nishitani, 2017-10-31

* file: /Users/bob/python/doing_math_with_python/linear_algebra.ipynb

Table of Contents¶

1 sympy

sympy¶

sympy tutorial
module references LUdecompositionの記述あり．

行列，ベクトルの生成¶

直接的に¶

In [1]:

from sympy import *

#init_session()
#init_printing(use_unicode=True)
init_printing()
list_a = [1,2]
list_b = [3,4]
A=Matrix([list_a, list_b])
A

Out[1]:

$\displaystyle \left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$

In [2]:

x,y = symbols('x,y')
v1=Matrix([x,y])
v1

Out[2]:

$\displaystyle \left[\begin{matrix}x\\y\end{matrix}\right]$

ゼロ行列，単位行列¶

In [3]:

zeros(2,3)

Out[3]:

$\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$

In [4]:

ones(2,3)

Out[4]:

$\displaystyle \left[\begin{matrix}1 & 1 & 1\\1 & 1 & 1\end{matrix}\right]$

In [5]:

eye(3)

Out[5]:

$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

In [6]:

diag(1,2,3)

Out[6]:

$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$

行列要素のとりだし，追加¶

In [7]:

A=Matrix([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
A[1,1]

Out[7]:

$\displaystyle 5$

In [8]:

A[:2,1:4]

Out[8]:

$\displaystyle \left[\begin{matrix}2 & 3\\5 & 6\end{matrix}\right]$

In [9]:

col_0 = A[:,0]
col_0

Out[9]:

$\displaystyle \left[\begin{matrix}1\\4\\7\end{matrix}\right]$

In [10]:

row_1 = A[1,:]
row_1

Out[10]:

$\displaystyle \left[\begin{matrix}4 & 5 & 6\end{matrix}\right]$

numpyではcolumn_stackを使って拡大係数行列を作った． sympyではrow_insert, col_insertを使う．

In [11]:

b = Matrix([1,2,3])
A.row_insert(3,b.T)

Out[11]:

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\\1 & 2 & 3\end{matrix}\right]$

In [12]:

A.col_insert(3,b)

Out[12]:

$\displaystyle \left[\begin{matrix}1 & 2 & 3 & 1\\4 & 5 & 6 & 2\\7 & 8 & 9 & 3\end{matrix}\right]$

行列・ベクトル間の単純な演算¶

次元の確認¶

In [13]:

A.shape

Out[13]:

$\displaystyle \left( 3, \ 3\right)$

In [14]:

print(v1.shape)
print(v1.T.shape)

(2, 1)
(1, 2)

転置(transpose)¶

In [15]:

A.T

Out[15]:

$\displaystyle \left[\begin{matrix}1 & 4 & 7\\2 & 5 & 8\\3 & 6 & 9\end{matrix}\right]$

In [16]:

v1.T

Out[16]:

$\displaystyle \left[\begin{matrix}x & y\end{matrix}\right]$

ドット積(内積)¶

In [17]:

A = Matrix([[1,2],[3,4]])
A*v1

Out[17]:

$\displaystyle \left[\begin{matrix}x + 2 y\\3 x + 4 y\end{matrix}\right]$

In [18]:

v1.T*v1 #v1*v1ではerror

Out[18]:

$\displaystyle \left[\begin{matrix}x^{2} + y^{2}\end{matrix}\right]$

In [19]:

v1.T*A

Out[19]:

$\displaystyle \left[\begin{matrix}x + 3 y & 2 x + 4 y\end{matrix}\right]$

In [20]:

v1.T*v1

Out[20]:

$\displaystyle \left[\begin{matrix}x^{2} + y^{2}\end{matrix}\right]$

In [21]:

v1*v1.T

Out[21]:

$\displaystyle \left[\begin{matrix}x^{2} & x y\\x y & y^{2}\end{matrix}\right]$

外積，outer, cross¶

outer(v1,v2)がSympy1.1.1では見つからない． crossは使える．

In [22]:

v1 = Matrix([1,1,3])
v2 = Matrix([1,2,-1])

In [23]:

v1.outer(v2)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-23-2fa213ccebb0> in <module>
----> 1 v1.outer(v2)

AttributeError: 'MutableDenseMatrix' object has no attribute 'outer'

In [24]:

v1.cross(v2)

Out[24]:

$\displaystyle \left[\begin{matrix}-7\\4\\1\end{matrix}\right]$

スカラー三重積¶

In [25]:

v3 = Matrix([-1,2,1])
v3.dot(v1.cross(v2))

Out[25]:

$\displaystyle 16$

行列操作¶

掃き出し, LU分解¶

In [ ]:

A=Matrix([[1,2,3],[4,5,6],[5,7,9]])
print(A.nullspace())
Aex = A.col_insert(3,b)
pprint(Aex)
Aex.rref()

rrefはreduced row echelon formを取り出してくれる．これはGaussの消去法で下んとこまで行って，さらに上んとこへ戻ってきたやつ．拡張行列もそのままやってくれるはず．もうすこし，試す必要があるが，一応MapleのLUDecompositionに近い結果を与えてくれる．

In [ ]:

A=Matrix([[1,2,3],[4,5,6],[7,9,8]])
L, U, P = A.LUdecomposition()
L, U, P

階数, rank¶

In [ ]:

print(A)

print(A.rank())

行列式, determinant¶

In [ ]:

c = Matrix([[1,1],[1,2]])
print(c.det())

a = Matrix([[1,2],[3,4]])
a.det()

課題¶

次の行列式の計算を確認せよ． $$ \left| \begin {array}{cccc} 1+a & 1 & 1 & 1\\ 1 & 1+b & 1 & 1\\ 1 & 1 & 1+c & 1\\ 1 & 1 & 1 & 1+d \end {array} \right| = abcd \left(1+\sum_{i=a}^d \frac{1}{i}\right) $$

In [ ]:

from sympy import *
a,b,c,d,x,y,z = symbols('a b c d x y z')

A = Matrix([[1+a,1,1,1],[1,1+b,1,1],[1,1,1+c,1],[1,1,1,1+d]]
)
pprint(A)
simplify(A.det())

逆行列¶

In [ ]:

A = Matrix([[1,2],[3,4]])
A.inv()

課題¶

次の連立方程式を解け． $$ \left\{ \begin {array}{cl} x+y+z&=0\\ ax+by+cz&=0\\ bcx+cay+abz&=(a-b)(b-c)(c-a) \end {array} \right. $$

In [ ]:

from pprint import pprint
from sympy import *
a,b,c,x,y,z = symbols('a b c x y z')

A = Matrix([[1,1,1],[a,b,c],[b*c,c*a,a*b]])
bb = Matrix([0,0,(a-b)*(b-c)*(c-a)])
print(A)
print(bb)
Ainv = A.inv()
res = Ainv * bb
pprint(simplify(res[0]))
pprint(simplify(res[1]))
pprint(simplify(res[2]))

固有値，固有ベクトル¶

sympyでは固有値・固有ベクトルは，eigenvects()で求められます．戻り値は，

[(固有値, 重複度, [固有ベクトル]), ...]

です．

In [ ]:

A = Matrix([[1,2],[2,1]])
A.eigenvects()

ベクトルの規格化¶

eigenvectsの固有ベクトルは規格化されていません．norm()を使うとかんたんに規格化してくれます．

In [ ]:

v1=Matrix(A.eigenvects()[0][2])
v2=Matrix(A.eigenvects()[1][2])

simplify(v1.T*v1)

In [ ]:

v_norm = v1.norm()
v_norm

In [ ]:

v1n = v1/v_norm
v1n

対角化¶

ちょいとわかりやすい行列を作成してみてみます．

In [ ]:

D = diag(1,2)
P = Matrix([[1,0],[1,1]]).T
A = P*D*P.inv()
A, P, D

In [ ]:

P,D = A.diagonalize()

P.inv()*A*P

In [ ]:

eigs = A.eigenvects()
v1 = eigs[0][-1][0]
v2 = eigs[1][-1][0]

Matrix([v1,v2]).reshape(2,2).T

対角化行列$P$をeigenvectsから作ろうとして苦労した．まず，固有ベクトルの取り出しが自動的にできない．重複しているときには，固有ベクトルは二つ出てくる．さらに，タプルからの取り出しがこれでいいいのか．．．さらに，ベクトルから行列を作るときに横置きのができない．うううんん．他のやり方をしたでやってみる．上では仕方がないので，一旦4x1のMatrixにしてそのあとreshape，さらにTransposeで．．．

In [ ]:

v1.col_insert(1,v2)

numpyにはravel()というフラット化の関数があるがsympyにはない．

結論：3つほど上でやっているdiagonalizeで取り出しましょう．重複があっても取り出してくれるんで．

課題¶

行列 $$ A\, := \, \left[ \begin {array}{ccc} 2&0&1\\ 0&3&0\\ 1&0&2 \end {array} \right] $$

を対角化する変換行列$P$を求め，対角化せよ．

In [ ]:

A = Matrix([[2,0,1],[0,3,0],[1,0,2]])
evs = A.eigenvects()
evs
#

In [ ]:

P,D = A.diagonalize()
P,D

In [ ]:

P.inv()*A*P

In [ ]: