2020年第1回のJDALのE資格を受験します。問題に応用数学が出るらしいので勉強した備忘録。以下①、②、③をまとめる予定です。
①逆行列計算
②特異値分解
③確率
※機械学習というか、統計モデル全般の考え方については下記記事が非常に参考になりました。
【訓練誤差と汎化誤差】学習・統計的推定は正しいのか？【過学習】
①逆行列計算
本項では以下の内容を記載します。
・2次正方行列の逆行列計算
・3次正方行列の逆行列計算
・numpyのlinalgを使った検算
任意行列の作成
今回例題に使用した行列は以下のコードで作成しました。
import numpy as np
n = 3 # 行列サイズ
A = np.random.randint(0,9,(n,n)) #0～9の整数でn×nの行列生成
2次正方行列の逆行列計算
A =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
A^{-1} = \frac{1}{det(A)}
\begin{pmatrix}
d & -b \\
-c & a
\end{pmatrix}
det(A) = ad-bc
◆例題
A =
\begin{pmatrix}
4 & 7 \\
6 & 9
\end{pmatrix}
A^{-1}=
\frac{1}{\rm 4・9-7・6}\
\begin{pmatrix}
9 & -7 \\
-6 & 4
\end{pmatrix}
=
\begin{pmatrix}
-3/2 & 7/6 \\
1 & -2/3
\end{pmatrix}
◆numpyのlinalgで検算
linalgはnumpyの線形代数モジュールです。詳しくは下記リンクを参照してください。
Linear algebra (numpy.linalg)
import numpy as np
A = np.array([[4,7],[6,9]])
# 逆行列の計算
invA = np.linalg.inv(A)
print(invA)
#[[-1.5 1.16666667]
# [1. -0.66666667]]
print(np.matmul(A,invA))
#[[ 1.00000000e+00 -3.33066907e-16]
# [ 0.00000000e+00 1.00000000e+00]]
# ちゃんと単位行列になってます
3次正方行列の逆行列計算
A =
\begin{pmatrix}
a_{11} & a_{12} & a_{13}\\
a_{21} & a_{22} & a_{23}\\
a_{31} & a_{32} & a_{33}
\end{pmatrix}
A^{-1} = \frac{1}{det(A)}
\begin{pmatrix}
Δ_{11} & Δ_{12} & Δ_{13}\\
Δ_{21} & Δ_{22} & Δ_{23}\\
Δ_{31} & Δ_{32} & Δ_{33}
\end{pmatrix}^{T}
det(A) = a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{21}a_{32} - a_{13}a_{22}a_{31} -a_{11}a_{23}a_{32}-a_{12}a_{21}a_{33}
Δ_{ij}：Aの i 行目と j 列目を除いた2×2行列の行列式を (−1)^{i+j}倍したもの(余因子)
◆行列式のイメージ
一般的にサラスの方法と呼ばれています。
◆余因子イメージ
◆例題
A =
\begin{pmatrix}
7 & 7 & 5\\
8 & 5 & 6\\
2 & 0 & 7
\end{pmatrix}
det(A)=
7・5・7+7・6・2+5・8・0-5・5・2-7・6・0-7・8・7 = -113
Δ_{11} = (-1)^{1+1}
\begin{vmatrix}
5 & 6 \\
0 & 7
\end{vmatrix}
=5・7-6・0=35
Δ_{12} = (-1)^{1+2}
\begin{vmatrix}
8 & 6 \\
2 & 7
\end{vmatrix}
=-8・7+6・2=-44\\
・・・以下略
A^{-1} =
\frac{-1}{113}\
\begin{pmatrix}
35 & -44 & -10\\
-49 & 39 & 14\\
17 & -2 & -21
\end{pmatrix}^{T}
=\frac{1}{113}\
\begin{pmatrix}
-35 & 49 & -17\\
44 & -39 & 2\\
10 & -14 & 21
\end{pmatrix}
◆numpyのlinalgで検算
import numpy as np
A = np.array([[7,7,5],[8,5,6],[2,0,7])
# 逆行列の計算
invA = np.linalg.inv(A)
print(invA)
#[[-0.30973451 0.43362832 -0.15044248]
# [ 0.38938053 -0.34513274 0.01769912]
# [ 0.08849558 -0.12389381 0.18584071]]
# 小数なのでわかりにくい
print(np.matmul(A,invA))
#[[ 1.00000000e+00 -1.80411242e-16 2.77555756e-17]
# [ 1.94289029e-16 1.00000000e+00 -5.55111512e-17]
# [ 4.16333634e-17 1.38777878e-17 1.00000000e+00]]
# ちゃんと単位行列になってはいる
逆行列が小数で答えが合ってるかわかりにくいため分数に直します。
まずpythonの組み込み関数のFractionを使って分数にしました。
from fractions import Fraction
print(Fraction(invA[0,0]))
# -2789840477132165/9007199254740992
丸め誤差の影響か元の分数に戻っていません。調べてもあまり小数⇒分数への変換コードが出てこないので、勉強も兼ねて自作しました。今回かなり力業で実装したので、もしもっと良い方法をご存じの方が居れば教えて頂きたいです。
以下処理の流れとソースです。
(1)あらかじめ分数のテーブルを作成
(2)逆行列要素との差が0に近いテーブル要素を抽出
(3)テーブル要素のindexから分数に再変換
# 分数のテーブルを用意
# table_size = 1000の時
# [[1/1,2/1,...,1000/1],[1/2,2/2,...,1000/2],...,[1/1000,2/1000,..,1000/1000]]
table_size = 3000
l_frac = np.array([[x/y for x in range(1,table_size+1)] for y in range(1,table_size+1)])
# 分数を浮動小数点に変換した際の誤差以下に設定
e = 1.0e-10
# 答えの分数を格納するndarray、後に"1/5"のように文字列で分数を入れるのでstrでキャスト
ans_flac = np.random.rand(len(A)**2).astype("str")
counter = 0
# 逆行列の要素一つずつの分数テーブルとの差を比較
for i in range(len(A)):
for j in range(len(A)):
if invA[i,j] != 0:
# 条件を満たす分数テーブルのindexを取得
# idx = (条件を満たす行のindex ,条件を満たす列のindex)
idx = np.where(abs(l_frac-abs(invA[i,j])) < e)
# 正負で条件分岐
if invA[i,j] < 0:
# 複数一致するのでindexの一番目(idx[][0])を選択しans_flacに格納
ans_flac[counter] = str(-(idx[1][0]+1)) + "/" + str(idx[0][0]+1)
counter += 1
else:
ans_flac[counter] = str(idx[1][0]+1) + "/" + str(idx[0][0]+1)
counter += 1
else:
# 逆行列要素が0の時
ans_flac[counter] = "0"
counter += 1
# 元の行列と同じ形状に変形
ans_flac = ans_flac.reshape([len(A),len(A)])
print(ans_flac)
#[['-35/113' '49/113' '-17/113']
# ['44/113' '-39/113' '2/113']
# ['10/113' '-14/113' '21/113']]
table_sizeから逸脱する分数が来るとエラーを吐くので、行列要素の桁が大きくなる場合はもっと大きいテーブルを用意する必要があります。
②特異値分解
この項では以下の内容を記載します。
・特異値分解の計算
・numpyのlinalgを使った検算
特異値の概念と使いどころは下記記事が参考になりました。
SVD（特異値分解）解説
上記記事内でも紹介されていますが、書籍ではゼロつくの自然言語処理編が詳しいです。
「ゼロから作る Deep Learning② 自然言語処理編」斎藤 康毅
詳細はリンク先に譲るとして特異値分解の定義と計算練習をします。
定義：m×nの行列Aについて、次の関係を満たすm×mの直行行列U,n×nの直行行列V、m×nの対角行列Sが存在する
A = USV^{T}
U,VはそれぞれAA^{T},A^{T}Aの固有ベクトルで作る行列で、ベクトルは左特異ベクトル、右特異ベクトルと呼ぶ。Sの対角成分は固有値の平方根であり、特異値と呼ばれる。
◆例題
A =
\begin{pmatrix}
1 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}
◆左側の行列U
まずは左特異ベクトルを求めて、行列Uを求めます。
AA^{T} =
\begin{pmatrix}
2 & 0\\
0 & 1
\end{pmatrix}\ 固有値をλと置くと\\
det(AA^{T}-λE) =
\begin{vmatrix}
2-λ & 0\\
0 & 1-λ
\end{vmatrix} = 0\\
∴λ = 1,2
λ = 2の時、固有ベクトルを(ux1,uy1)と置くと
\begin{pmatrix}
2-λ & 0\\
0 & 1-λ
\end{pmatrix}
\begin{pmatrix}
u_{x1}\\
u_{y1}
\end{pmatrix} =
\begin{pmatrix}
0 & 0\\
0 & -1
\end{pmatrix}
\begin{pmatrix}
u_{x1}\\
u_{y1}
\end{pmatrix}
= 0\\
\\
\left\{
\begin{array}{ll}
u_{x1}・0 + u_{y1}・0 = 0 \\
u_{x1}・0 - u_{y1}・1 = 0
\end{array}
\right.
これが常に成立する時、
∴u_{y1} = 0
したがってλ=2の時の固有ベクトルは、ベクトルの長さを1に規格化すると、次のようになります。
\begin{pmatrix}
u_{x1}\\
u_{y1}
\end{pmatrix}=
\begin{pmatrix}
1\\
0
\end{pmatrix}・・・(1)
同様にλ=1の時、固有ベクトルを(ux2,uy2)と置くと
\begin{pmatrix}
2-λ & 0\\
0 & 1-λ
\end{pmatrix}
\begin{pmatrix}
u_{x2}\\
u_{y2}
\end{pmatrix} =
\begin{pmatrix}
1 & 0\\
0 & 0
\end{pmatrix}
\begin{pmatrix}
u_{x2}\\
u_{y2}
\end{pmatrix}
= 0\\
\\
\left\{
\begin{array}{ll}
u_{x2}・1 + u_{y2}・0 = 0 \\
u_{x2}・0 + u_{y2}・0 = 0
\end{array}
\right.\\
∴u_{y1} = 0
固有ベクトルは下記です。
\begin{pmatrix}
u_{x2}\\
u_{y2}
\end{pmatrix}=
\begin{pmatrix}
0\\
1
\end{pmatrix}・・・(2)
(1)、(2)からUが求められます。
U =
\begin{pmatrix}
u_{x1} & u_{x2}\\
u_{y1} & u_{y2}
\end{pmatrix}=
\begin{pmatrix}
1 & 0\\
0 & 1
\end{pmatrix}
◆右側の行列V
次に右特異ベクトルを求めて、行列Vを求めます。
A^{T}A =
\begin{pmatrix}
1 & 1 & 0\\
1 & 1 & 0\\
0 & 0 & 1
\end{pmatrix}\ 固有値をλと置くと\\
det(AA^{T}-λE) =
\begin{vmatrix}
1-λ & 1 & 0\\
1 & 1-λ & 0\\
0 & 0 & 1-λ
\end{vmatrix} = 0\\
(1-λ)(λ-2)λ = 0\\
∴λ = 0,1,2
途中計算は省きましたが、サラスの方法で行列式を計算できます。
λ = 2の時、固有ベクトルを(vx1,vy1,vz1)と置くと
\begin{pmatrix}
1-λ & 1 & 0\\
1 & 1-λ & 0\\
0 & 0 & 1-λ
\end{pmatrix}
\begin{pmatrix}
v_{x1}\\
v_{y1}\\
v_{z1}
\end{pmatrix} =
\begin{pmatrix}
-1 & 1 & 0\\
1 & -1 & 0\\
0 & 0 & -1
\end{pmatrix}
\begin{pmatrix}
v_{x1}\\
v_{y1}\\
v_{z1}
\end{pmatrix}
= 0\\
\\
\left\{
\begin{array}{ll}
-1・v_{x1} + 1・v_{y1} + 0・v_{z1} = 0 \\
1・v_{x1} - 1・v_{y1} + 0・v_{z1} = 0 \\
0・v_{x1} + 0・v_{y1} - 1・v_{z1} = 0
\end{array}
\right.
これが常に成立する時、
∴v_{x1} = v_{y1},v_{z1} = 0
したがって固有ベクトルは
\begin{pmatrix}
v_{x1}\\
v_{y1}\\
v_{z1}
\end{pmatrix}=
\begin{pmatrix}
1/\sqrt{2}\\
1/\sqrt{2}\\
0
\end{pmatrix}・・・(3)
λ=1,0の時も同様に計算できます。ここでは答えだけ記載します。
λ=1の時
\begin{pmatrix}
v_{x2}\\
v_{y2}\\
v_{z2}
\end{pmatrix}=
\begin{pmatrix}
0\\
0\\
1
\end{pmatrix}・・・(4)\\
λ=0の時
\begin{pmatrix}
v_{x3}\\
v_{y3}\\
v_{z3}
\end{pmatrix}=
\begin{pmatrix}
1/\sqrt{2}\\
-1/\sqrt{2}\\
0
\end{pmatrix}・・・(5)
(3),(4),(5)よりVが求められます。
V =
\begin{pmatrix}
v_{x1} & v_{x2} & v_{x3}\\
v_{y1} & v_{y2} & v_{y3}\\
v_{z1} & v_{z2} & v_{z3}
\end{pmatrix}=
\begin{pmatrix}
1/\sqrt{2} & 0 & 1/\sqrt{2}\\
1/\sqrt{2} & 0 & -1/\sqrt{2}\\
0 & 1 & 0
\end{pmatrix}
◆特異値行列S
固有値の平方根を大きい順に対角に並べます。
S =
\begin{pmatrix}
\sqrt{2} & 0 & 0\\
0 & 1 & 0
\end{pmatrix}
※三つ目の固有値0の扱いについては勉強不足で理解しきれていません。わかり次第更新します。
◆検算
ここでもnumpyのlinalgモジュールを使って答え合わせをします。linalg.svd()という関数でU,S,Vを取得できます。
import numpy as np
A = np.array([[1,1,0],[0,0,1]])
# Vは転置が計算されます
U,S,V_T = np.linalg.svd(A, full_matrices = True)
S = np.diag(S) #特異値は[λ1,λ2]のように一次元なので対角化
S = np.insert(S,2,[0,0],axis = 0) #次元合わせ
print("U=\n",U)
print("S=\n",S)
print("V=\n",V.T) #転置の転置で元のVを表示
print("A=\n",np.dot(np.dot(U,S),V))
# U=
# [[1. 0.]
# [0. 1.]]
# S=
# [[1.41421356 0. 0. ]
# [0. 1. 0. ]]
# V=
# [[ 0.70710678 0. 0.70710678]
# [ 0.70710678 0. -0.70710678]
# [ 0. 1. 0. ]]
# A=
# [[1. 1. 0.]
# [0. 0. 1.]]
# 元の行列と一致
その他線形代数の基礎は下記記事も参考になりました。
機械学習の基礎２ 線形代数メモ
③確率
この項では以下の内容を記載します。
・条件付き確率
・ベイズの定理
・ベルヌーイ分布
・二項分布
・マルチヌーイ(カテゴリカル)分布
・多項分布
確率分布に関しては下記サイトが非常にわかりやすく、参考にさせていただきました。
ややこしい離散分布に関するまとめ
条件付き確率
AとB二人が10本中3本あたりの入ったくじを一回ずつ引く場合を考えます。くじは引いた後元に戻すことにします。
あたり確率を$P(A)$、はずれ確率を$P(\overline{A})$、のように表記しベン図と表にまとめると以下のようになります。
この時Aがあたりを引き、Bもあたりを引く条件付き確率を$P(B \mid A)$と置くと、
P(B \mid A) =\frac{P(A∩B)}{P(A)}・・・(6)
となります。下記図のイメージです。
※AとBが同時に成り立つ同時確率$P(A∩B)$は下記図のイメージです。
ベイズの定理
条件付き確率の式(6)により次の関係も成り立ちます。
P(A \mid B) =\frac{P(A∩B)}{P(B)}・・・(7)
(6)、(7)から$P(A∩B)$を消去した式がベイズの定理です。
P(B \mid A) =\frac{P(A \mid B) P(B)}{P(A)}・・・(8)
使いどころは別途更新します。
ベルヌーイ分布
結果が2通りしかない事象を考えます。良く例に出されるのはコインを投げて裏が出るか1度だけ試す場合です。表の出る確率を$μ$とすると、裏表の確率は以下にまとめられます。
\left\{
\begin{array}{ll}
表の出る確率 μ = 0.5 \\
裏の出る確率 1-μ = 0.5 \\
\end{array}
\right.
またここでコイン表の事象を$X = 1$、裏の事象を$X = 0$と置くと、コイン投げの確率分布$Bern(X \mid μ)$は
Bern(X \mid μ) =
μ^{X}(1-μ)^{1-X} \\(X=0,1)
となります。グラフで示すとこんな感じです。
# ベルヌーイ分布の描画
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
mu = 0.5
x = np.array([0,1])
y = (mu**x)*(1-mu)**(1-x)
plt.scatter(x,y)
また、期待値$E(X) = μ$、分散$σ(X) = μ(1-μ)$です。
二項分布
コイン投げを複数回繰り返す場合の確率分布です。N回コインを投げて、m回表が出る二項分布$Bin(m \mid\ μ,N)$は
Bin(m \mid\ μ,N)=\frac{N!}{m!(N-m)!}μ^{m}(1-μ)^{N-m}
となります。
また、期待値$E(X) = Nμ$、分散$σ(X) = Nμ(1-μ)$です。
マルチヌーイ(カテゴリカル)分布
ベルヌーイ分布の拡張です。結果がn通りある事象を考えます。良く例に出されるのはサイコロを振って出目を試す場合です。各目の出る確率を$μ_{k}$とすると、マルチヌーイ分布$Cat(X\mid\textbf{μ})$は以下になります。
Cat(X\mid\textbf{μ})=\prod_{k=1}^{n}\mu_k^{X_{k}}
多項分布
マルチヌーイ分布の拡張です。結果がn通りある事象をN回施行する場合の確率分布です。多項分布$Multi(\textbf{m}\mid\textbf{μ},N)$は以下になります。
Multi(\textbf{m}\mid\textbf{μ},N)= \frac{N!}{m_{1}!...m_{k}!}\prod_{k=1}^{n}\mu_k^{m_{k}}
正規分布
一般的に一番よく使われている分布。
f(x) = \frac{1}{\sqrt{2\pi\sigma}}\exp{\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}}
$\mu$は平均、$\sigma^{2}$は分散。
$Z = (x-\mu)/\sigma$と置くと、標準正規分布になります。
f(x) = \frac{1}{\sqrt{2\pi}}\exp{\left\{-\frac{Z^2}{2}\right\}}
おまけ
ノルム
$L_{p}$ノルムの定義
{\| {\bf x} \|_p = (\ |x_1|^p + |x_2|^p + \cdots + |x_n|^p )^{1/p}}
正則化の文脈で出てきます。$L_{2}$ノルムはユークリッド距離、$L_{1}$ノルムはマンハッタン距離とも呼ばれます。過学習を防ぐweight decayとして損失関数に付与します。
【機械学習】LPノルムってなんだっけ？
[DL]weight decayって何？