線形代数(Linear algebra)

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

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

sympy

行列,ベクトルの生成

直接的に

In [2]:
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[2]:
$$\left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$$
In [9]:
x,y = symbols('x,y')
v1=Matrix([x,y])
v1
Out[9]:
$$\left[\begin{matrix}x\\y\end{matrix}\right]$$

ゼロ行列,単位行列

In [10]:
zeros(2,3)
Out[10]:
$$\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$
In [11]:
ones(2,3)
Out[11]:
$$\left[\begin{matrix}1 & 1 & 1\\1 & 1 & 1\end{matrix}\right]$$
In [12]:
eye(3)
Out[12]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$
In [13]:
diag(1,2,3)
Out[13]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 2 & 0\\0 & 0 & 3\end{matrix}\right]$$

行列要素のとりだし,追加

In [14]:
A=Matrix([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
A[1,1]
Out[14]:
$$5$$
In [15]:
A[:2,1:4]
Out[15]:
$$\left[\begin{matrix}2 & 3\\5 & 6\end{matrix}\right]$$
In [16]:
col_0 = A[:,0]
col_0
Out[16]:
$$\left[\begin{matrix}1\\4\\7\end{matrix}\right]$$
In [17]:
row_1 = A[1,:]
row_1
Out[17]:
$$\left[\begin{matrix}4 & 5 & 6\end{matrix}\right]$$

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

In [18]:
b = Matrix([1,2,3])
A.row_insert(3,b.T)
Out[18]:
$$\left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\\7 & 8 & 9\\1 & 2 & 3\end{matrix}\right]$$
In [19]:
A.col_insert(3,b)
Out[19]:
$$\left[\begin{matrix}1 & 2 & 3 & 1\\4 & 5 & 6 & 2\\7 & 8 & 9 & 3\end{matrix}\right]$$

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

次元の確認

In [46]:
A.shape
Out[46]:
$$\left ( 3, \quad 3\right )$$
In [49]:
print(v1.shape)
print(v1.T.shape)
(3, 1)
(1, 3)

転置(transpose)

In [20]:
A.T
Out[20]:
$$\left[\begin{matrix}1 & 4 & 7\\2 & 5 & 8\\3 & 6 & 9\end{matrix}\right]$$
In [21]:
v1.T
Out[21]:
$$\left[\begin{matrix}x & y\end{matrix}\right]$$

ドット積(内積)

In [23]:
A = Matrix([[1,2],[3,4]])
A*v1
Out[23]:
$$\left[\begin{matrix}x + 2 y\\3 x + 4 y\end{matrix}\right]$$
In [29]:
v1.T*v1 #v1*v1ではerror
Out[29]:
$$\left[\begin{matrix}x^{2} + y^{2}\end{matrix}\right]$$
In [24]:
v1.T*A
Out[24]:
$$\left[\begin{matrix}x + 3 y & 2 x + 4 y\end{matrix}\right]$$
In [25]:
v1.T*v1
Out[25]:
$$\left[\begin{matrix}x^{2} + y^{2}\end{matrix}\right]$$
In [26]:
v1*v1.T
Out[26]:
$$\left[\begin{matrix}x^{2} & x y\\x y & y^{2}\end{matrix}\right]$$

外積,outer, cross

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

In [30]:
v1 = Matrix([1,1,3])
v2 = Matrix([1,2,-1])
In [31]:
v1.outer(v2)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-31-419c7c24b3ff> in <module>()
----> 1 v1.outer(v2)

/Users/bob/anaconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py in __getattr__(self, attr)
   3214         else:
   3215             raise AttributeError(
-> 3216                 "%s has no attribute %s." % (self.__class__.__name__, attr))
   3217 
   3218     def integrate(self, *args):

AttributeError: MutableDenseMatrix has no attribute outer.
In [32]:
v1.cross(v2)
Out[32]:
$$\left[\begin{matrix}-7\\4\\1\end{matrix}\right]$$

スカラー三重積

In [33]:
v3 = Matrix([-1,2,1])
v3.dot(v1.cross(v2))
Out[33]:
$$16$$

行列操作

掃き出し, LU分解

In [42]:
A=Matrix([[1,2,3],[4,5,6],[5,7,9]])
print(A.nullspace())
Aex = A.col_insert(3,b)
pprint(Aex)
Aex.rref()
[Matrix([
[ 1],
[-2],
[ 1]])]
⎡1  2  3  1⎤
⎢          ⎥
⎢4  5  6  2⎥
⎢          ⎥
⎣5  7  9  3⎦
Out[42]:
$$\left ( \left[\begin{matrix}1 & 0 & -1 & - \frac{1}{3}\\0 & 1 & 2 & \frac{2}{3}\\0 & 0 & 0 & 0\end{matrix}\right], \quad \left [ 0, \quad 1\right ]\right )$$

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

In [43]:
A=Matrix([[1,2,3],[4,5,6],[7,9,8]])
L, U, P = A.LUdecomposition()
L, U, P
Out[43]:
$$\left ( \left[\begin{matrix}1 & 0 & 0\\4 & 1 & 0\\7 & \frac{5}{3} & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 2 & 3\\0 & -3 & -6\\0 & 0 & -3\end{matrix}\right], \quad \left [ \right ]\right )$$

階数, rank

In [26]:
print(A)

print(A.rank())
Matrix([[1, 2, 3], [4, 5, 6], [7, 9, 8]])
3

行列式, determinant

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

a = Matrix([[1,2],[3,4]])
a.det()
1
Out[27]:
$$-2$$

課題

次の行列式の計算を確認せよ. $$ \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 [59]:
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())
⎡a + 1    1      1      1  ⎤
⎢                          ⎥
⎢  1    b + 1    1      1  ⎥
⎢                          ⎥
⎢  1      1    c + 1    1  ⎥
⎢                          ⎥
⎣  1      1      1    d + 1⎦
Out[59]:
$$a b c d + a b c + a b d + a c d + b c d$$

逆行列

In [60]:
A = Matrix([[1,2],[3,4]])
A.inv()
Out[60]:
$$\left[\begin{matrix}-2 & 1\\\frac{3}{2} & - \frac{1}{2}\end{matrix}\right]$$

課題

次の連立方程式を解け. $$ \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 [28]:
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]))
Matrix([[1, 1, 1], [a, b, c], [b*c, a*c, a*b]])
Matrix([[0], [0], [(-a + c)*(a - b)*(b - c)]])
-b + c
a - c
-a + b

固有値,固有ベクトル

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

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

です.

In [13]:
A = Matrix([[1,2],[2,1]])
A.eigenvects()
Out[13]:
$$\left [ \left ( -1, \quad 1, \quad \left [ \left[\begin{matrix}-1\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 1, \quad \left [ \left[\begin{matrix}1\\1\end{matrix}\right]\right ]\right )\right ]$$

ベクトルの規格化

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

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

simplify(v1.T*v1)
Out[14]:
$$\left[\begin{matrix}2\end{matrix}\right]$$
In [15]:
v_norm = v1.norm()
v_norm
Out[15]:
$$\sqrt{2}$$
In [18]:
v1n = v1/v_norm
v1n
Out[18]:
$$\left[\begin{matrix}- \frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\end{matrix}\right]$$

対角化

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

In [19]:
D = diag(1,2)
P = Matrix([[1,0],[1,1]]).T
A = P*D*P.inv()
A, P, D
Out[19]:
$$\left ( \left[\begin{matrix}1 & 1\\0 & 2\end{matrix}\right], \quad \left[\begin{matrix}1 & 1\\0 & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 0\\0 & 2\end{matrix}\right]\right )$$
In [20]:
P,D = A.diagonalize()

P.inv()*A*P
Out[20]:
$$\left[\begin{matrix}1 & 0\\0 & 2\end{matrix}\right]$$
In [21]:
eigs = A.eigenvects()
v1 = eigs[0][-1][0]
v2 = eigs[1][-1][0]

Matrix([v1,v2]).reshape(2,2).T
Out[21]:
$$\left[\begin{matrix}1 & 1\\0 & 1\end{matrix}\right]$$

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

In [22]:
v1.col_insert(1,v2)
Out[22]:
$$\left[\begin{matrix}1 & 1\\0 & 1\end{matrix}\right]$$

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

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

課題

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

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

In [23]:
A = Matrix([[2,0,1],[0,3,0],[1,0,2]])
evs = A.eigenvects()
evs
#
Out[23]:
$$\left [ \left ( 1, \quad 1, \quad \left [ \left[\begin{matrix}-1\\0\\1\end{matrix}\right]\right ]\right ), \quad \left ( 3, \quad 2, \quad \left [ \left[\begin{matrix}0\\1\\0\end{matrix}\right], \quad \left[\begin{matrix}1\\0\\1\end{matrix}\right]\right ]\right )\right ]$$
In [24]:
P,D = A.diagonalize()
P,D
Out[24]:
$$\left ( \left[\begin{matrix}-1 & 0 & 1\\0 & 1 & 0\\1 & 0 & 1\end{matrix}\right], \quad \left[\begin{matrix}1 & 0 & 0\\0 & 3 & 0\\0 & 0 & 3\end{matrix}\right]\right )$$
In [25]:
P.inv()*A*P
Out[25]:
$$\left[\begin{matrix}1 & 0 & 0\\0 & 3 & 0\\0 & 0 & 3\end{matrix}\right]$$
In [ ]: