Numpy & Scipy - 1.4 行列の基本操作 (2)

hstack / vstack

• np.hstack((tuple)), np.vstack((tuple))：2D array, 1D array、二つが混合された形式で組み合わせが可能です。
• 両方とも deep copy です。

• (1) 2D array の場合

• 2x3 2x2 hstack

a = np.array([[1,2,3], [4,5,6]], dtype=np.float64)
b = np.array([[-1,-2], [-3,-4]], dtype=np.float64) 

new_mat = np.hstack( (a,b) ) # 入力時 tuple 形式であることに注意

new_mat = np.hstack( (a,b,b,b) ) # 複数組み合わせ可能

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

# 水平に行列が組み合わされたことを確認できます。
1.0, 2.0, 3.0,-1.0,-2.0
4.0, 5.0, 6.0,-3.0,-4.0

# 2x5 matrix
(2, 5)

• 3x3 2x2 hstack

a = np.array([[1,2,3],[4,5,6]], dtype=np.float64)
b = np.array([[-1,-2,-3],[-4,-5,-6],[-7,-8,-9]], dtype=np.float64)

new_mat = np.hstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

# 2x3 3x3 は 2row と 3row hstack 組み合わせが不可能
ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 2 and the array at index 1 has size 3

• 2x3 2x2 vstack

a = np.array([[1,2,3], [4,5,6]], dtype=np.float64)
b = np.array([[-1,-2], [-3,-4]], dtype=np.float64) 

new_mat = np.vstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

# 2x3 2x2 は vstack 組み合わせができないため次のようなエラーが発生します。

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 3 and the array at index 1 has size 2

• 2x3 3x3 vstack

a = np.array([[1,2,3],[4,5,6]], dtype=np.float64)
b = np.array([[-1,-2,-3],[-4,-5,-6],[-7,-8,-9]], dtype=np.float64)

new_mat = np.vstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

1.0, 2.0, 3.0
4.0, 5.0, 6.0
-1.0,-2.0,-3.0
-4.0,-5.0,-6.0
-7.0,-8.0,-9.0

(5, 3)

(2) 1D array の場合

• hstack
• 1D array ベクトルがより長くなります。

a = np.array([1,2,3], dtype=np.float64)
b = np.array([4,5,6], dtype=np.float64)

new_mat = np.hstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

1.0, 2.0, 3.0, 4.0, 5.0, 6.0

(6,)

• vstack
• 行列が作られます。

a = np.array([1,2,3], dtype=np.float64)
b = np.array([4,5,6], dtype=np.float64)

new_mat = np.vstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

1.0, 2.0, 3.0
4.0, 5.0, 6.0

(2, 3)

• (3) 2D array と 1D array 混合された形式

• hstack
• hstack は無条件にエラーが発生します。混合された形式は vstack のみ可能です。

a = np.array([[1,2,3],[4,5,6]],dtype=np.float64)
b = np.array([7,8,9], dtype=np.float64)

new_mat = np.hstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

• vstack
• 2x3 2D array と 1x3 1D array を組み合わせ

a = np.array([[1,2,3],[4,5,6]],dtype=np.float64)
b = np.array([7,8,9], dtype=np.float64)

new_mat = np.vstack( (a,b) ) # 入力時 tuple 形式であることに注意

prt(new_mat, fmt="%0.1f", delimiter=",")
print()
print(new_mat.shape)

1.0, 2.0, 3.0
4.0, 5.0, 6.0
7.0, 8.0, 9.0

(3, 3)

transpose method

• 文字通り行列を transpose 転置する関数

• swallow copy です。

• (1) transpose property

c = np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=np.float64)

d = c.T

prt(c, fmt="%0.1f", delimiter=",")
print()
prt(d, fmt="%0.1f", delimiter=",")

# c matrix
1.0, 2.0, 3.0
4.0, 5.0, 6.0
7.0, 8.0, 9.0

# d matrix (transpose of c)
1.0, 4.0, 7.0
2.0, 5.0, 8.0
3.0, 6.0, 9.0

• swallow copy 例示

• (2) transpose method

• 1D array を transpose() するとそのまま 1D array で出ます。

e = np.array([1,2,3], dtype=np.float64)

f = e.transpose()

g = np.copy(e.transpose()) # np.copy を通じて deep copy が可能

prt(e, fmt="%0.1f", delimiter=",")
print(e.shape)
print()
prt(f, fmt="%0.1f", delimiter=",")
print(f.shape)

1.0, 2.0, 3.0
(3,)

1.0, 2.0, 3.0
(3,)

real / imag / conjugate

• real, imag は各行列 entry たちの実数部分、虚数部分を持ってきて返します。
• imag の場合、虚数部分を持ってきて実数行列 or ベクトルを返します。

• 両方とも swallow copy です。

• conjugate は deep copy です。

a = np.array([[1-2j, 3+1j, 1], [1+2j, 2-1j, 7]], dtype=np.complex128)

a_real = r.real

# prt(a, fmt="%0.1f", delimiter=",")
# print()
# prt(a_real, fmt="%0.1f", delimiter=",")

a_imag = r.imag # 虚数部分を実数で返します

# prt(a, fmt="%0.1f", delimiter=",")
# print()
# prt(a_imag, fmt="%0.1f", delimiter=",")


a_conj = r.conjugate() # conjugate を取ると虚数部分に - が掛けられて出ます
# deep copy

# prt(a, fmt="%0.1f", delimiter=",")
# print()
# prt(a_conj, fmt="%0.1f", delimiter=","

# a matrix
( 1.0-2.0j),( 3.0+1.0j),( 1.0+0.0j)
( 1.0+2.0j),( 2.0-1.0j),( 7.0+0.0j)

# real
1.0, 3.0, 1.0
1.0, 2.0, 7.0

# imag
-2.0, 1.0, 0.0
 2.0,-1.0, 0.0

# conjugate
( 1.0+2.0j),( 3.0-1.0j),( 1.0+0.0j)
( 1.0-2.0j),( 2.0+1.0j),( 7.0+0.0j)

Multiplication

• Scalar Multiplication
• 順序が関係ありません result = r * A = A * r

A = np.array([[1,2,1],[2,1,3],[1,3,1]], dtype=np.float64)
scalar = 5.0

result = scalar * A

prt(A, fmt="%0.1f", delimiter=",")
print()
prt(result, fmt="%0.1f", delimiter=",")

1.0, 2.0, 1.0
2.0, 1.0, 3.0
1.0, 3.0, 1.0

5.0, 10.0, 5.0
10.0, 5.0, 15.0
5.0, 15.0, 5.0

• Matrix Multiplication

• 順序が非常に重要です。数学的行列乗算と同じです。

• @ operator を使用する方法と
• np.matmul() 関数が存在します。

A = np.array([[1,2,3], [3,2,1]], dtype=np.float64)
B = np.array([[2,1],[1,2],[-3, 3]], dtype=np.float64)

result = A @ B
result2 = np.matmul(A, B)

prt(A, fmt="%0.1f", delimiter=",")
print()
prt(B, fmt="%0.1f", delimiter=",")
print()
prt(result, fmt="%0.1f", delimiter=",")
print()
prt(result2, fmt="%0.1f", delimiter=",")

1.0, 2.0, 3.0
3.0, 2.0, 1.0

2.0, 1.0
1.0, 2.0
-3.0, 3.0

# A @ B
-5.0, 14.0
 5.0, 10.0

# np.matmul(A,B)
-5.0, 14.0
 5.0, 10.0

• Matrix-Vector Multiplication
• 同様に result = A @ u , result = np.matmul(A,u) を使用して計算が可能です。
• dot product ではありません！

A = np.array([[1,2,1],[2,1,3],[1,3,1]], dtype=np.float64)
u = np.array([5,1,3], dtype=np.float64)

Au = np.matmul(A, u)

prt(A, fmt="%0.1f", delimiter=",")
print()
prt(u, fmt="%0.1f", delimiter=",")
print()
prt(Au, fmt="%0.1f", delimiter=",")

1.0, 2.0, 1.0
2.0, 1.0, 3.0
1.0, 3.0, 1.0

5.0, 1.0, 3.0

10.0, 20.0, 11.0

Inner Product : vdot

• result = np.vdot(u, v)

• ここで u, v はベクトルです。

• Inner Product を行うと次のように linear combination 形式で表現されます。

real vector : $\quad \mathbf{u} \cdot \mathbf{v} = u_1v_1 + u_2v_2 + \dots + u_nv_n$
complex vector : $ \quad \mathbf{u} \cdot \mathbf{v} = u_1\bar{v}_1 + u_2\bar{v}_2 + \dots + u_n\bar{v}_n = \bar{u_1v_1} + \bar{u_2v_2} + \dots + \bar{u_nv_n}$

• 例示..

u = np.array([1,1,1,1], dtype=np.float64)
v = np.array([-1,1,-1,1], dtype=np.float64)

vdot = np.vdot(u,v)

prt(u, fmt="%0.1f", delimiter=",")
print()
prt(v, fmt="%0.1f", delimiter=",")
print()
print(vdot)

 1.0, 1.0, 1.0, 1.0

-1.0, 1.0,-1.0, 1.0

0.0

Norm

• 私たちが以前まで線形代数学で学んだベクトルの norm は \(\lVert \mathbf{x} Vert _2 = \left( \sum_{i=1}^n x_i^2 ight)^{1/2} \quad l_2 \mbox{-vector norm (Euclidean)}\) このように定義し、これは $l_2$ norm と呼びます。これ以外にも多くのベクトル norm が存在します。

vector norm

\[\lVert \mathbf{x} Vert _1 = \sum_{i=1}^n \lvert x_i vert \quad l_1 \mbox{-vector norm}\] \[\lVert \mathbf{x} Vert _2 = (\sum_{i=1}^n \lvert x_i^2 vert)^{1/2} \quad l_2 \mbox{-vector norm (Euclidean)}\] \[\lVert \mathbf{x} Vert _{\infty} = \max_{1 \le i \le n} \lvert x_i vert\]

• そして同様に行列の norm も存在しますが、ベクトル norm のように $l_1, l_2, l_{\infty}$ norm が存在します。

matrix norm

\[\lVert A Vert _1 = \max_{1 \le j \le n} \sum_{i=1}^m \lvert a_{ij} vert \quad l_1 \mbox{-matrix norm}\] \[\lVert A Vert _2 = \sigma _{\max} \le \left( \sum_{i=1}^m \sum_{j=1}^n \lvert a_{ij} vert ^2 ight)^{1/2} \quad l_2 \mbox{-matrix norm (spectral)}\] \[\lVert A Vert _{\infty} = \max_{1 \le i \le m} \sum_{j=1}^n \lvert a_{ij} vert \quad l_{\infty} \mbox{-matrix norm}\]

• numpy でこの linalg.norm を使用して求めるには scipy ライブラリを使用する必要があります。

from scipy import linalg

• 例示

w = np.array([1,1], dtype=np.float64)

norm = linalg.norm(w, 2)

prt(w, fmt="%0.1f", delimiter=",")
print()
print(norm)

 1.0, 1.0

 1.4142135623730951

• また $l_{\infty}$ norm を求めるには inf をパラメータとして入れれば良いです。

norm = linalg.norm(a, 2, inf)

slicing - 行列またはベクトルの一部だけを抜き出す

• swallow copy です。

• 例示

A = np.array([[1,2,3,4,5], [6,7,8,9,10], [11,12,13,14,15], [-1,-2,-3,-4,-5]])
sub_A = A[1:3, 1:4]

• row : 1~2, column : 1~3 を持ってくるという意味です。
• つまり : の後ろにあるインデックス -1 まで持ってくるという意味です！
• 1:3 を入力すると実際には 1:2 まで持ってきます。

sub_A = A[0:5, 2:6]

# 以下と同一です

sub_A = A[ : , 2: ]

• : は最初から最後まで
• 2: は 2 の後ろから全部

• row, column 側に slicing operator : を付けないと 1D array で返します。
• np.reshape を通じて 2D array に変換できます。

sub_A1 = A[1:2, 1:4] # 2D array で返却
sub_A2 = A[1, 1:4] # 1D array で返却

sub_A3 = A[1: , 1:2] #2D array で返却
sub_A4 = A[1:, 1] #1D araay で返却

sub_A5 = np.reshape(sub_A3, (3,1)) # 1D array -> 2D array

#sub_A1
7.0, 8.0, 9.0
(1, 3)

#sub_A2
7.0, 8.0, 9.0
(3,)

#sub_A3
 7.0
 12.0
-2.0
(3, 1)

#sub_A4
7.0, 12.0,-2.0
(3,)

#sub_A5
 7.0
 12.0
-2.0
(3, 1)

Example 1.

$A = \begin{bmatrix} 1 & 2 \ 3 & 4 \end{bmatrix} \quad \mathbf{x} = \begin{bmatrix} 5 \ 6 \end{bmatrix} $
$A, \mathbf{x}$ をそれぞれ変数に入れて 2D, 1D array で保存し、1D array は transpose ができないことを勘案して quadratic form である $\mathbf{x}^TA\mathbf{x}$ を求めてみる。
1) 2D array に変換してみる (np.reshape)
2) np.vdot 使用してみる

import numpy as np
from print_lecture import print_custom as prt
from scipy import linalg

A = np.array([[1,2],[3,4]], dtype=np.float64)
x = np.array([5,6], dtype=np.float64)

# 1) 
x_t = np.reshape(x, (2,1))

result1 = np.matmul(np.matmul(x_t.T, A), x)

print(result1)

# 2) np.vdot 

result2 = np.vdot(x, A @ x)

print(result2)

# 1番 結果
319.0

# 2番 結果
319.0