e資格対策(1) Numpyで線形代数

はじめ

インフラネットワークにもデータサイエンスを取り入れて活用できるんじゃないのかな? というのが始まりで、まずは”AIを知る”意味で昨年末にg検定認定を受けた。 次は”実装できる”を目標にe資格取得の勉強を始めたので、その備忘録。

インフラネットワークとコラボできるかは、、、うーん未知。でも時代とともにインフラだってインフラエンジニアに求められるものだって変わっていくはず、いや、変えていきたいよね。

一回目 Numpyで線形代数

前提準備として、PCにanacondaのインストールが済んでいること。以下からDL可能。 今回の操作はすべてWindowsコマンドプロンプトのみ。

www.python.jp

やること

  • Numpyのバージョン確認
  • ランダムな数値で行列を作成する
  • 逆行列を求める
  • 行列式の値を求める
  • ベクトルの内積を求める
  • ベクトルの外積を求める
  • 単位行列を作成する
  • 全ての成分が0の行列を作成する
  • ベクトルのL1、L2のノルムを求める
  • 行列の積を求める
  • 固有値を求める
  • 固有ベクトルを求める
  • 転置行列を求める
  • 直行行列の定義を確認する
  • 特異値分解(SVD)・・・(2)に続く
用語

逆行列・・・行列の逆数。特徴として(行列)(逆行列)=単位行列.行列Aに対してA-1で表す
単位行列・・・(identity matrix)1が斜め右下に伸びる対角行列.Iで表す
転置(T)・・・行列の縦横逆転した行列.行列Aに対してATで表す
正方行列・・・2×2や3×3など、正方形になる要素の行列
正則行列・・・正方行列の中で逆行列が1つ存在する行列(1つのみが普通)
固有値分解・・・正方行列を分解すること.正方行列Aに対してA=VΛV-1
特異値分解・・・非正方行列を分解すること.非正方行列MにおいてM=USV-1
左特異ベクトル・・・単位ベクトルの直行行列.非正方行列MにおいてMMT固有値分解して求める.Uで表す
右特異クトル・・・単位ベクトルの直行行列.非正方行列MにおいてMMT固有値分解して求める.Vで表す
固有値・・・ (英: eigenvalue) 正方行列でλで表す.
特異値・・・σで表す
固有ベクトル・・・ (英: eigenvector) 正方行列で固有値が並んだ対角行列ベクトルV
単位ベクトル・・・大きさが1のベクトル
直行行列・・・(orthogonal matrix)転置行列と逆行列が等しくなる正方行列.Rで表す
対称行列・・・対角成分で2分すると対称となる行列。対称行列は直行行列で対角化可能
対角化・・・行列から対角行列を取り出す.固有値分解のΛは対角行列なので要するに固有値分解のこと。対角行列=P-1APを書き換えるとA=P対角行列P-1でA=VΛV-1と同じ
対角行列・・・左列の上から対角線上に数値が並び他は0の行列

Numpyのバージョン確認

※以後、不明な点はバージョンに合ったマニュアルを参照する。

NumPy 1.17.4 Release Notes — NumPy v1.19.dev0 Manual

>>>
>>> import numpy as np
>>>
>>> print(np.__version__)
1.17.4
>>>
ランダムな数値で行列を作成する

numpy.random.randintで求めることができる。0から9までの整数値を使って3×3の行列を作成する例。

>>> ma4=np.random.randint(0, 9, (3, 3))
>>>
>>> ma4
array([[7, 7, 2],
       [5, 3, 0],
       [2, 7, 7]])
>>>
>>> ma4=np.random.randint(0, 9, (3, 3))
>>>
>>> inv(ma4)
array([[ 0.        ,  0.33333333, -1.66666667],
       [ 0.5       ,  0.        , -1.5       ],
       [-0.5       ,  0.        ,  2.5       ]])
>>>
逆行列を求める

numpy.linalg.inv()で求めることができる。

>>> from numpy.linalg import inv
>>>
>>> ma1=np.array([[1,1],[2,1]])
>>>
>>> ma1
array([[1, 1],
       [2, 1]])
>>>
>>> inv(ma1)
array([[-1.,  1.],
       [ 2., -1.]])
>>>

※ad-bd=0のためエラーとなる例

>>> ma2=np.array([[1,1],[1,1]])
>>>
>>> ma2
array([[1, 1],
       [1, 1]])
>>>
>>> inv(ma2)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<__array_function__ internals>", line 5, in inv
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 551, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 97, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix
>>>


※正方行列でないためエラーとなる例

>>> ma3= np.array([[2,2,2,2],[1,-1,1,-1],[-1,1,-1,1]])
>>>
>>> ma3
array([[ 2,  2,  2,  2],
       [ 1, -1,  1, -1],
       [-1,  1, -1,  1]])
>>>
>>> inv(ma3)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<__array_function__ internals>", line 5, in inv
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 546, in inv
    _assertNdSquareness(a)
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 213, in _assertNdSquareness
    raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square
>>>
>>>
行列式の値を求める

行列式(英: determinant)で求めることができる。正方行列の要素同士の演算で求められるスカラー(ベクトルではない)。‘numpy.linalg.det‘で計算できる。 例は3×3の行列式の値と、2×3の行列ではエラーとなる例。

>>> a = np.array([[1,2,1],[4,3,1],[0,1,0]])
>>>
>>> a
array([[1, 2, 1],
       [4, 3, 1],
       [0, 1, 0]])
>>>
>>> np.linalg.det(a)
3.0000000000000004
>>>

>>> a = np.array([[1,2],[4,3],[0,1]])
>>>
>>>
>>> a
array([[1, 2],
       [4, 3],
       [0, 1]])
>>>
>>> np.linalg.det(a)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<__array_function__ internals>", line 5, in det
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 2122, in det
    _assertNdSquareness(a)
  File "C:\Users\miwapinn\AppData\Local\Programs\Python\Python38-32\lib\site-packages\numpy\linalg\linalg.py", line 213, in _assertNdSquareness
    raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square
ベクトルの内積を求める

※numpy.inner(a,b)で求めることができる。

>>> a=np.array([1,2,1])
>>> b=np.array([3,1,1])
>>>
>>> np.inner(a,b)
6
ベクトルの外積を求める

※numpy.outer(a,b)で求めることができる。

>>> b=np.array([3,1,1])
>>> a=np.array([1,2,1])
>>>
>>> np.outer(a,b)
array([[3, 1, 1],
       [6, 2, 2],
       [3, 1, 1]])
>>>
単位行列を作成する

※対角成分が1でそれ以外が全て0の行列、numpy.identity(次数)で求めることができる。3×3と5×5の例

>>>
>>> np.identity(3)
array([[1., 0., 0.]
       [0., 1., 0.],
       [0., 0., 1.]])
>>>
>>> a=np.identity(5)
>>>
>>> a
array([[1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0.],
       [0., 0., 1., 0., 0.],
       [0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 1.]])
>>>
全ての成分が0の行列を作成する

np.zeros()で作成できる。3×1と4×4の例を記載

>>> np.zeros(3)
array([0., 0., 0.])
>>>
>>> np.zeros([4,4])
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])
>>>
ベクトルのL1、L2のノルムを求める

※ノルムはLA.norm()で計算できる。以下のLA.norm(a, 1)の1はL1ノルム。2の場合はL2ノルムで数字を指定しない場合のデフォルト。L2ノルムは各要素の2乗和の平均となる。

f:id:pinn38:20200509172251p:plain
抜粋 https://www.python.jp/install/anaconda/index.html

>>> from numpy import linalg as LA
>>>
>>> a=np.array([[0,0],[2,2]])
>>>
>>> a
array([[0, 0],
       [2, 2]])
>>>
>>> LA.norm(a, 1)
2.0
>>>
>>> LA.norm(a, 2)
2.8284271247461903
>>>
>>> LA.norm(a)
2.8284271247461903
>>>
>>>
行列の積を求める

※いくつか方法がある。1つはnp.matrixで変数を行列として定義すれば*で積は求められる。または、np.matmulnp.dotで計算する。または@で計算する。 以下のa,bはnp.matrixで定義し、c,dはnp.arrayで定義し各計算をしてみた例

>>> a = np.array([[1,2],[2,3]])
>>> b = np.array([[0,1],[1,1]])
>>>
>>> c = np.matrix([[1,2],[2,3]])
>>> d = np.matrix([[0,1],[1,1]])
>>>
>>> a * b
array([[0, 2],
       [2, 3]])
>>>
>>> c * d
matrix([[2, 3],
        [3, 5]])
>>>
>>> a @ b
array([[2, 3],
       [3, 5]])
>>>
>>> c @ d
matrix([[2, 3],
        [3, 5]])
>>>
>>> np.matmul(a,b)
array([[2, 3],
       [3, 5]])
>>>
>>> np.matmul(c,d)
matrix([[2, 3],
        [3, 5]])
>>>
>>> np.dot(a,b)
array([[2, 3],
       [3, 5]])
>>>
>>> np.dot(c,d)
matrix([[2, 3],
        [3, 5]])
>>>
固有値を求める

※numpy.linalg パッケージをLAとしてimportしLA.eig()で求める。
(例)

>>> a = np.array([[2,0],[1,6]])
>>>
>>> LA.eigvals(a)
array([6., 2.])
>>>
固有ベクトルを求める(固有値と同時に求める)

(例)
w・・・固有値
v・・・固有ベクトル
固有ベクトルは大きさ1に正則化して出力される(固有ベクトルは必ずしも正則化しなくてもOK)

>>> a = np.array([[2,0],[1,6]])
>>>
>>> w,v = LA.eig(a)
>>>
>>> w
array([6., 2.])
>>>
>>> v
array([[ 0.        ,  0.9701425 ],
       [ 1.        , -0.24253563]])
>>>
転置行列を求める

※行と列の要素を入れかえた行列、aの転置g行列をbとするとb = a.Tで計算できる

>>> a=np.array([[3,2,0],[0,2,0],[0,0,1]])
>>>
>>> a
array([[3, 2, 0],
       [0, 2, 0],
       [0, 0, 1]])
>>>
>>> b = a.T
>>>
>>> b
array([[3, 0, 0],
       [2, 2, 0],
       [0, 0, 1]])
>>>
直行行列の定義を確認する

対称行列は直交行列により対角化可能<<<★ポイント
※直行行列Rの定義
RT=R-1
RRT=RTR=I
積: (RS)T (RS) = (RS)(RS)T = I
行列式|R| = +- 1
固有値λ = +- 1

ここまでのnumpyの関数を使って確認する。Rは直行行列、Iは単位行列を表す。

まず、対称行列Aを定義

>>> A = np.array([[2,1,0],[1,1,1],[0,1,2]])
>>>
>>> A
array([[2, 1, 0],
       [1, 1, 1],
       [0, 1, 2]])

以下計算によりAは固有値Sで対角化可能.
ここで、Aが対称行列の場合(今回はそう)固有値ベクトルRは直行行列であり、 R-1=RTが確認できる.
※見にくいので小数点以下6桁で丸めている

>>> A = np.array([[2,1,0],[1,1,1],[0,1,2]])
>>>
>>> f,R = LA.eig(A)
>>>
>>> np.round(R.T,6)
array([[ 0.408248, -0.816497,  0.408248],
       [ 0.707107,  0.      , -0.707107],
       [ 0.57735 ,  0.57735 ,  0.57735 ]])
>>>
>>> np.round(inv(R),6)
array([[ 0.408248, -0.816497,  0.408248],
       [ 0.707107,  0.      , -0.707107],
       [ 0.57735 ,  0.57735 ,  0.57735 ]])
>>>

次にRRT = RTR = Iを以下で確認.
※0にマイナス符号がついているのは6桁で丸めているため.

>>> np.round(R@RT,6)
array([[ 1.,  0.,  0.],
       [ 0.,  1., -0.],
       [ 0., -0.,  1.]])
>>>
>>> np.round(RT@R,6)
array([[ 1., -0., -0.],
       [-0.,  1.,  0.],
       [-0.,  0.,  1.]])
>>>

行列式を計算し、|R|=+-1であることを確認.

>>> A = np.array([[2,1,0],[1,1,1],[0,1,2]])
>>>
>>> f,R = LA.eig(A)
>>>
>>> np.abs(np.linalg.det(R))
1.0
>>>

Rの固有値の絶対値が1であることを確認.
λ = +- 1

>>> A = np.array([[2,1,0],[1,1,1],[0,1,2]])
>>>
>>> f,R = LA.eig(A)
>>>
>>> LA.eigvals(R)
array([-0.00720072+0.99997407j, -0.00720072-0.99997407j,
        1.        +0.j        ])
>>>
>>> np.abs(LA.eigvals(R))
array([1., 1., 1.])
>>>