線形代数(Linear algebra) scipy版

cc by Shigeto R. Nishitani, 2017-2018

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

pythonでの線形代数のlibraryはいくつもあるように見える. 混乱しがちなんで,それをいくつか分類してみた. 間違ってたら教えて.

  • numpy.linalg
  • scipy.linalg
  • sympy

scipyはnumpyのwrapperっぽい.NumpyとScipyによると,

  1. NumPy ⊂ SciPy ということのようだ
  2. numpy で提供されている機能はそのまま, scipy でも提供されている
  3. なので scipy だけで押し通しても良さそうだが, 世の中の説明は numpy が主流なので, それに合わせて, 基本は numpy, scipy だけで 提供されている機能は scipy を使う

と明言されている.これが一番混乱がなさそう.

sympyは代数計算向きなんで,表記がだいぶ違う.別章として提供.

出力精度の制御

行列の出力はとても醜くなる場合がある. ほぼ0なのにとか,次元が増えてとかで.そのときこいつが役立つ.

In [1]:
import numpy as np
from pprint import pprint
import scipy.linalg as linalg

a = np.array([[2,5], [4,1]])
l,P = np.linalg.eig(a)
pprint(l)
pprint(P)

np.set_printoptions(precision=3, suppress=True)

pprint(l)
pprint(P)
array([ 6., -3.])
array([[ 0.78086881, -0.70710678],
       [ 0.62469505,  0.70710678]])
array([ 6., -3.])
array([[ 0.781, -0.707],
       [ 0.625,  0.707]])

行列,ベクトルの生成

In [2]:
import numpy as np
import scipy as sp
list_a = [1,2]
list_b = [3,4]
np_a = np.array([list_a, list_b])
print(np_a)

sp_a = sp.array([list_a, list_b])
print(sp_a)

np_m = np.matrix([list_a, list_b])
print(np_m)
[[1 2]
 [3 4]]
[[1 2]
 [3 4]]
[[1 2]
 [3 4]]
In [3]:
list_a = [1,2]

np_a_1 = np.array(list_a)
print(np_a_1)

np_v = np.array([list_a])
print(np_v)
# np.vectorはない.
[1 2]
[[1 2]]

ゼロ行列,単位行列

In [10]:
np.zeros((3,3))
Out[10]:
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
In [8]:
np.diag([1,2,3])
Out[8]:
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
In [9]:
np.identity(3)
Out[9]:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

転置(transpose)

In [7]:
print(np_a.transpose())
print(sp_a.transpose())
print(np_m.transpose())
[[1 3]
 [2 4]]
[[1 3]
 [2 4]]
[[1 3]
 [2 4]]
In [11]:
np_v.transpose()
Out[11]:
array([[1],
       [2]])

行列,ベクトルの演算

ドット積

$$ \left[ \begin {array}{cc} 1&2\\ 3&4 \end {array} \right] \left[ \begin {array}{c} 1\\ 2 \end {array} \right] $$

を計算したいときは,

In [4]:
np.dot(np_a,np.transpose(np_v))
Out[4]:
array([[ 5],
       [11]])
In [5]:
np_v.dot(np_a)
Out[5]:
array([[ 7, 10]])

これだと, $$ \left[ \begin {array}{cc} 1&2 \end {array} \right] \left[ \begin {array}{cc} 1&2\\ 3&4 \end {array} \right] $$

を計算することになる.

ベクトルの距離 $$ \left[ \begin {array}{cc} 1&2 \end {array} \right] \left[ \begin {array}{cc} 1\\ 2 \end {array} \right] $$ は

In [6]:
np_v.dot(np_v)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-491742187a7f> in <module>()
----> 1 np_v.dot(np_v)

ValueError: shapes (1,2) and (1,2) not aligned: 2 (dim 1) != 1 (dim 0)

ではダメで,

In [7]:
np_v.dot(np_v.transpose())
Out[7]:
array([[5]])
In [8]:
np_v.transpose().dot(np_v)
Out[8]:
array([[1, 2],
       [2, 4]])

と順番を間違うと悲惨.単に長さをだけなら$||v||$の代わりにnp.linalg.norm(v)も使える.

In [9]:
np.linalg.norm(np_v,2)**2
Out[9]:
5.0000000000000009

ただし,normなんでそのほかの指数も指定できる.

In [10]:
np.linalg.norm(np_v,1)**2
Out[10]:
4.0

外積,outer, cross

In [30]:
v1 = np.array([1,1,3])
v2 = np.array([1,2,-1])
In [31]:
np.outer(v1,v2)
Out[31]:
array([[ 1,  2, -1],
       [ 1,  2, -1],
       [ 3,  6, -3]])
In [32]:
np.cross(v1,v2)
Out[32]:
array([-7,  4,  1])

スカラー三重積

In [14]:
v3 = np.array([-1,2,1])
np.dot(np.cross(v1,v2),v3)
Out[14]:
16

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

In [33]:
a=np.array([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
a[1,1]
Out[33]:
5
In [34]:
a[:2,1:4]
Out[34]:
array([[2, 3],
       [5, 6]])
In [35]:
col_0 = a[:,0]
print(col_0)
[1 4 7]
In [36]:
row_1 = a[1,:]
print(row_1)
[4 5 6]
In [37]:
b = np.array([1,2,3])
np.vstack((a,b))
Out[37]:
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9],
       [1, 2, 3]])

hstackはa,bが同じ次元でないとダメとのお達し.拡大係数行列を手軽に作れるわけではなさそう.column_stackを使うとできた.

In [38]:
np.column_stack((a,b))
Out[38]:
array([[1, 2, 3, 1],
       [4, 5, 6, 2],
       [7, 8, 9, 3]])

掃き出し, LU分解

In [43]:
import scipy.linalg

a = np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])

P, L, U = scipy.linalg.lu(a)

print(P)
print(L)
print(U)
[[1. 0. 0.]
 [0. 0. 1.]
 [0. 1. 0.]]
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 0. 1.]]
[[ 1. -1. -1.]
 [ 0.  2.  0.]
 [ 0.  0.  2.]]
In [45]:
b = np.array([1,-1,-1])
ab = np.column_stack((a,b))
P, L, U = scipy.linalg.lu(ab)
print(P)
print(L)
print(U)
[[1. 0. 0.]
 [0. 0. 1.]
 [0. 1. 0.]]
[[1. 0. 0.]
 [1. 1. 0.]
 [1. 0. 1.]]
[[ 1. -1. -1.  1.]
 [ 0.  2.  0. -2.]
 [ 0.  0.  2. -2.]]

階数

In [42]:
print(a)

print(np.linalg.matrix_rank(a))
print(np.rank(a)) #deprecatedなんでやめろって,
# 関数名とかlibの区分けに統一性がないよね...
[[1 2 3]
 [4 5 6]
 [7 8 9]]
2
2
/Users/syasin/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  after removing the cwd from sys.path.

逆行列

In [9]:
a = np.array([[1,2],[3,4]])
scipy.linalg.inv(a)
Out[9]:
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

行列式

In [8]:
c = np.array([[1,1],[1,2]])
print(np.linalg.det(c))

a = np.array([[1,2],[3,4]])
scipy.linalg.det(a)
1.0
Out[8]:
-2.0

連立方程式の解

逆行列を用いて,連立方程式の解を求めるには次の通り. 例えば,連立方程式が次のような場合, $$ x - y -z = 1 \\ x - y +z = -1 \\ x + y -z = -1 $$

In [10]:
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
scipy.linalg.inv(a).dot(b)
Out[10]:
array([-1., -1., -1.])

課題

次の連立方程式を解け. $$ \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 [26]:
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

固有値,固有ベクトル

In [11]:
a = np.array([[1,1],[0,2]])
l, P = scipy.linalg.eig(a)
print(l)
print(P)
[ 1.+0.j  2.+0.j]
[[ 1.          0.70710678]
 [ 0.          0.70710678]]

ベクトルの規格化

In [12]:
# eigの固有ベクトルは規格化されている
p_norm = np.linalg.norm(P[:,1])
p_norm
# 規格化はこいつをつかって...
Out[12]:
0.99999999999999989

対角化

Pに固有ベクトルが入っているので,対角化は$P^{-1} A P$で

In [13]:
np.dot(np.dot(np.linalg.inv(P),a),P)
#np.dot(np.dot(np.linalg.inv(P),a),P)
Out[13]:
array([[ 1.,  0.],
       [ 0.,  2.]])

課題

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

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

In [30]:
import numpy as np

A = np.array([[2,0,1],[0,3,0],[1,0,2]])
l, P = np.linalg.eig(A)
print(l)
print(P)

np.dot(np.dot(np.linalg.inv(P),A),P)
[ 3.  1.  3.]
[[ 0.70710678 -0.70710678  0.        ]
 [ 0.          0.          1.        ]
 [ 0.70710678  0.70710678  0.        ]]
Out[30]:
array([[  3.00000000e+00,   4.80650139e-17,   0.00000000e+00],
       [ -1.01465364e-17,   1.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   3.00000000e+00]])
In [31]:
import numpy as np
from math import pi

np.linspace(0,2*pi,9)
Out[31]:
array([ 0.        ,  0.78539816,  1.57079633,  2.35619449,  3.14159265,
        3.92699082,  4.71238898,  5.49778714,  6.28318531])
In [ ]: