pythonでの線形代数のlibraryはいくつもあるように見える. 混乱しがちなんで,それをいくつか分類してみた. 間違ってたら教えて.
scipyはnumpyのwrapperっぽい.NumpyとScipyによると,
と明言されている.これが一番混乱がなさそう.
sympyは代数計算向きなんで,表記がだいぶ違う.別章として提供.
行列の出力はとても醜くなる場合がある. ほぼ0なのにとか,次元が増えてとかで.そのときこいつが役立つ.
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]])
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]]
<ipython-input-3-6b5cbc232797>:8: DeprecationWarning: scipy.array is deprecated and will be removed in SciPy 2.0.0, use numpy.array instead sp_a = sp.array([list_a, list_b])
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]]
print(np_v.shape)
print(np_m.shape)
(1, 2) (2, 2)
np.zeros((3,3))
array([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]])
np.diag([1,2,3])
array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
np.identity(3)
array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]])
print(np_a.transpose())
print(sp_a.transpose())
print(np_m.transpose())
[[1 3] [2 4]] [[1 3] [2 4]] [[1 3] [2 4]]
np_v.transpose()
array([[1], [2]])
print(np_v.shape)
print(np_v.transpose().shape)
(1, 2) (2, 1)
を計算したいときは,
np.dot(np_a,np.transpose(np_v))
array([[ 5], [11]])
np_v.dot(np_a)
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] $$ は
np_v.dot(np_v)
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-13-2f795847ce59> in <module> ----> 1 np_v.dot(np_v) ValueError: shapes (1,2) and (1,2) not aligned: 2 (dim 1) != 1 (dim 0)
ではダメで,
np_v.dot(np_v.transpose())
array([[5]])
np_v.transpose().dot(np_v)
array([[1, 2], [2, 4]])
と順番を間違うと悲惨.単に長さだけなら$||v||$の代わりにnp.linalg.norm(v)も使える.
np.linalg.norm(np_v)**2
5.000000000000001
ただし,normなんでそのほかの指数も指定できる.
np.linalg.norm(np_v,1)**2
4.0
v1 = np.array([1,1,3])
v2 = np.array([1,2,-1])
np.outer(v1,v2)
array([[ 1, 2, -1], [ 1, 2, -1], [ 3, 6, -3]])
np.cross(v1,v2)
array([-7, 4, 1])
v3 = np.array([-1,2,1])
np.dot(np.cross(v1,v2),v3)
16
a=np.array([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
a[1,1]
5
a[:2,1:4]
array([[2, 3], [5, 6]])
col_0 = a[:,0]
print(col_0)
[1 4 7]
row_1 = a[1,:]
print(row_1)
[4 5 6]
b = np.array([1,2,3])
np.vstack((a,b))
array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [1, 2, 3]])
hstackはa,bが同じ次元でないとダメとのお達し.拡大係数行列を手軽に作れるわけではなさそう.column_stackを使うとできた.
np.column_stack((a,b))
array([[1, 2, 3, 1], [4, 5, 6, 2], [7, 8, 9, 3]])
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.]]
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.]]
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.
a = np.array([[1,2],[3,4]])
scipy.linalg.inv(a)
array([[-2. , 1. ], [ 1.5, -0.5]])
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
-2.0
逆行列を用いて,連立方程式の解を求めるには次の通り. 例えば,連立方程式が次のような場合, $$ x - y -z = 1 \\ x - y +z = -1 \\ x + y -z = -1 $$
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
scipy.linalg.inv(a).dot(b)
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. $$
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
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]]
# eigの固有ベクトルは規格化されている
p_norm = np.linalg.norm(P[:,1])
p_norm
# 規格化はこいつをつかって...
0.99999999999999989
Pに固有ベクトルが入っているので,対角化は$P^{-1} A P$で
np.dot(np.dot(np.linalg.inv(P),a),P)
#np.dot(np.dot(np.linalg.inv(P),a),P)
array([[ 1., 0.], [ 0., 2.]])
行列 $$ A\, = \, \left[ \begin {array}{ccc} 2&0&1\\ 0&3&0\\ 1&0&2 \end {array} \right] $$
を対角化する変換行列$P$を求め,対角化せよ.
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. ]]
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]])
import numpy as np
from math import pi
np.linspace(0,2*pi,9)
array([ 0. , 0.78539816, 1.57079633, 2.35619449, 3.14159265, 3.92699082, 4.71238898, 5.49778714, 6.28318531])