データサイエンスに必須な高校数学・大学数学の基礎知識を一瞬で復習しつつ特異値分解を理解する【行列・ベクトル・線型代数】

こんにちは。

最近このブログでシリーズ化しているペンギンデータセットを使ったデータサイエンス入門ですが、次に扱おうとしている「主成分分析」の記事の前に、一度基本的な数学の知識を復習しておかないと、進むのが難しくなってきました。

というわけで、今回はデータサイエンスや統計学で特に重要な線型代数・行列・ベクトルの要点を、公式と概要を見ながらさらっと復習してしまおう、という記事です。というか自分が忘れてるので振り返りたいだけ。

特に、「主成分分析」に必要な「特異値分解」の理解に辿り着くまでに必要な知識を抑えられるような構成で組み立てています。特異値分解というのはこういうやつです。

うん、まだ何か全然わからない。けれどこの記事を読めばわかるようになるはず。

皆さんも「こんなのやったな〜」みたいな軽いノリで眺めるだけでも、良い復習になると思います。ぜひノスタルジー(?)に浸りつつ、今後のデータサイエンスに活用していきましょう。ビッグデータと言えど、所詮は巨大な行列やベクトルなので、基礎が大事。

理系の方で、一度は習ったことある方向けに書いているので、本当に本当の基本的なところはすっ飛ばして書きます。ご容赦ください。

プログラミングに活かせるように、本文中にはPythonのコードも併せて書いておきます。ライブラリにはNumPyを使います。

最初に、インポートして$2 \times 2$の行列(2次元配列)の$A$を作っておきます。(本文はこの時点で何言っているかわかるレベルの人向けに書いてます。)

import numpy as np
A = np.array([[2,3],[4, -1]])

また、NumPyやSciPyの英文ドキュメントの理解に役立つように、できるだけ英語名も記載していきます。これ大事。そのままメソッドの名前にも使われてたりするので。

それでは見ていきましょう。

単位行列 identity matrix

対角成分が全て1の行列。$E$とか$I$とかで使う。Einheitsというドイツ語派か、Identityという英語派か、で記号が分かれる。

どちらから掛けても元の行列$A$と同じになる特性を持つ。

I = np.eye(3) # 3行3列の単位行列
# array([[1., 0., 0.],
#        [0., 1., 0.],
#        [0., 0., 1.]])

参考 : NumPy | numpy.eye

行列の積(ドット積、点乗積) dot product

次の連立方程式は、

行列で表すとこうなる。

というわけで、積$AB$を計算したい場合、$A$の列数と$B$の行数が一致していないと計算できない。$A$が2列なら、$B$は2行でないと計算できない。

ちなみに、$A^{2}$といったようにべき乗を計算したいときは、$A$の列数と行数が同じ、つまり正方行列でなくてはならない。

NumPyでは積の計算はdotを使う。もしくはmatrixにしてから*でもよい。この結果のように$AB$と$BA$は必ずしも一致しない。

A = np.array([[2,3],[4, -1]])
B = np.array([1,2])
np.dot(A,B) # array([8, 2])
np.dot(B,A) # array([10,  1])

べき乗するときはmatrixにしないとできない。arrayのままでは各要素で演算されるため結果が異なる。

np.matrix(A)**2
# matrix([[16,  3],
#         [ 4, 13]])
A**2
# array([[ 4,  9],
#        [16,  1]])

また、Python3.5以降であれば@演算子が使える。知ってると書きやすいし読みやすいので超便利。ただ、知らないと何コレってなる。

A@B # array([8, 2])

行列式 determinant

行列ではなくスカラーになります。

行列の積の行列式は、行列式の積に一致する。

また、転置行列$A^{T}$の行列式は、もとの行列の行列式と同じになる。

np.linalg.det(A) # -14.000000000000004
np.linalg.det(A.T) # .Tメソッドで転置できる。detの結果は同じ-14。

ちなみにlinalgはLinear algebraの略で、線型代数のこと。

参考 : NumPy | Linear algebra (numpy.linalg)

正方行列 square matrix

$2\times2$や$100\times100$のように、行数と列数が一致する行列のこと。

似たような言葉に正則行列(regular matrix)もある。こちらは、行列式が0ではない場合の行列のことで、つまり逆行列を持つ正方行列のこと。可逆行列(invertible matrix)とか非特異行列(non-singular matrix)とか呼ばれたりもする。

逆行列 inverse matrix

正方行列$A$について、

が成り立つとき、$B = A^{-1}$と書ける。これを$A$の逆行列という。

2行2列の場合は、次の公式で求まる。

公式を一般化すると余因子(cofactor)とかいうマニアックな知識が必要になるので、ここではひとまず割愛。

逆行列は、Pythonで簡単に求められます。そう、NumPyならね。

np.linalg.inv(A)
# array([[ 0.07142857,  0.21428571],
#        [ 0.28571429, -0.14285714]])
np.matrix(A).I # matrixはIメソッドが使える。結果はmatrixで返る。

ベクトルの内積 inner product

次の関係で表されるスカラーです。内積が0なら直交しているという性質は超重要。

ちなみにベクトルは太い小文字のアルファベットで書くのが通例ですが、わかりにくいので本記事では高校数学の慣習に倣って$\vec{a}$のように矢印を乗せて書きます。

また、 $\vec{a} = \begin{pmatrix}a_1\\a_2\end{pmatrix}$、 $\vec{b} = \begin{pmatrix}b_1\\b_2\end{pmatrix}$ の場合、内積の計算はこうなる。

同じベクトルなら角度が一致しているので$cos$が1になり、$|\vec{a}|^{2}$。

Pythonで計算するとこう。ベクトル同士の内積は、ドット積と同義なので、.dot()を使います。ということは、もちろん@演算子も使えます。

a = np.array([1,2])
b = np.array([3,4])
np.dot(a,b) # 11
np.dot(a,a) # 5
a@b # 11
a@a # 5

基底 basis と成分 component と次元 dimension

ちょっと小難しい用語の話ですが、NumPyとかのドキュメントに出てくるので一応おさらい。

基底とは、線型空間(linear space)$V$の中で、一次独立かつ$V$の要素を網羅的に記述できるベクトルの組。

たとえば、$xy$座標系だと、$ax+by$で座標内のどの点も記述できるので、このとき$\vec{x}$と$\vec{y}$が基底。基底の長さがそれぞれ1でかつ直交していたら正規直交基底(orthonornal basis)とかも言ったりします。

空間内のある点$a\vec{x}+b\vec{y}$を表すときは、$(a,b)$と書き、これを基底の成分(component)と言う。

ちなみに、一次独立(linearly independent)というのは一次結合$a\vec{x}+b\vec{y}$したときに、零ベクトルにするためには0倍しないといけない、つまり$(a,b)=(0,0)$しかないということ。逆に0以外の組み合わせもあれば、一次従属という。懐かしい。

また、$V$の中で基底となるベクトルの数を次元(dimension)と言います。先の$xy$座標なら、基底の数は$x$と$y$の2つなので、2次元。まあ、普通のことをカッコつけて言ってるだけです。数学的には$dimV = 2$と書くらしい。

次元をPythonで書くとこう。配列の次元が返ってきます。

A.ndim #2

階数、ランク rank

行列の階数(ランク)とは、行列$A$の中で、一次独立である行(または列)の数のこと。

2次の正方行列なら最大ランク2で、$n$次の正方行列の場合は最大ランク$n$となる。

このような行列の場合、2行目は1行目の2倍なので、一次独立な行は1。つまりランクは1。

np.linalg.matrix_rank(A) # 2
np.linalg.matrix_rank(np.array([[1,2],[2,4]])) # 1

線型写像 linear mapping

線型変換(linear transformation)ともいう。しゃぞーってなんすか？のアレ。

線型空間$V$で、ある行列やベクトルで変換された後の座標だったりベクトルのことで、写像を$f$とすると、以下の性質を持つ。

たとえば$xy$座標の平面空間で、ある座標$(1,2)$が行列$A$によって変換されると、座標は$(10,1)$に移動する。このように$xy$の任意の点が移動していった先が、写像。

ちなみに$\begin{pmatrix}1\\2\end{pmatrix}$による$x$座標の写像は$y=2x$の直線になるといったように、１次関数もある意味線型写像だったりする。これが一番わかり易い。(と思う)

固有値 eigenvalue と固有ベクトル eigenvector

ぶっちゃけここからの章を書きたくて、この記事を作り始めました。

正方行列$A$について、次の式を満たす$\vec{x}$が固有ベクトルで、$\lambda$で表されるスカラーが固有値。

$A$が正方行列の場合、$\vec{x}$を線型変換、つまり同じ空間内で別のベクトルに変換していることになります。それが、スカラー倍で表されるということは、方向は変わらないけれど長さが変わる変換が行われているということです。このような特殊なベクトルを固有ベクトルといいます。

ということは、固有値は正方行列でないと求められないのか！と思いがちですが、後ほどこれを正方行列ではない行列に一般化した「特異値」というものも出てきます。そっちが本番。

固有方程式(特性方程式) eigen equation

固有ベクトルを知りたい！というときは、まず以下の方程式で固有値を求めます。固有値を求めるこの方程式を、固有方程式といいます。

これを満たす$\lambda$すべてが、$A$の固有値です。そして、以下の式を満たす$\vec{x}$が、固有ベクトル。

NumPyで求めるには、linalg.eig()メソッド。

w, v = np.linalg.eig(A)
w # array([ 4.27491722, -3.27491722])　これがAの固有値で、2つある
v # これがAの固有ベクトルで、2つの固有値に対してそれぞれ存在する。
# array([[ 0.79681209, -0.49436913],
#       [ 0.60422718,  0.86925207]])

ちなみにvは正規化されて長さが1になっています。

たしかにwとvは上の固有方程式や関係式を満たしていそうなことが、以下からわかります。

I = np.eye(2)
np.linalg.det(A - w[0]*I) # 0.0 : 0になった！
np.linalg.det(A - w[1]*I) # 0.0

np.dot(A - w[0]*I, v[:,0]) # array([0., 0.]) : 零ベクトルになった！
np.dot(A - w[1]*I, v[:,1]) # array([0., 0.])

固有ベクトルvは、1個目が横ではなく縦に並んでいるので、v[:,0]というスライスを使うか、v.T[0]のように転置して取り出すのに注意。

参考 : NumPy | numpy.linalg.eig

対角和 trace

ここはちょっと余談ですが、さらに固有値についてマニアックな性質になると、正方行列$A$の対角の和は、固有値の総和と一致します。

対角というのは単位行列でいう1が並ぶところです。

np.trace(A) # Aの対角和 = 1
w[0] + w[1] = 1.0

対称行列 symmetric matrix

自身の転置と一致する正方行列。対角以外の成分が対称になっている。

こんなやつ。

対角行列 diagonal matrix

こちらは対称行列と名前がややこしいけれど、対角以外の成分がすべて0の行列。

こんなやつ。

こちらはとても重要な性質で、固有値は対角成分と一致します。つまり、この行列の場合、固有値は2と-1。

NumPyだと、.diag()で対角成分を抜き出した対角行列を作れます。上の例は実は何度も登場している行列$A$君の対角成分。

np.diag(A) # array([ 2, -1])
np.diag(np.diag(A))
# array([[ 2,  0],
#       [ 0, -1]])
np.linalg.eig(np.diag(np.diag(A)))
# (array([ 2., -1.]),
# array([[1., 0.],
#        [0., 1.]]))

たしかに固有値が対角成分と一致している。

あれ？つまり対角行列が作れれば、固有値求めるのってめっちゃ楽なのでは？というのが続く話。

固有値分解 eigen value decomposition

固有値分解とは、行列$A$を次の形に変形すること。真ん中に、固有値$\lambda$を対角成分に持つ対角行列$\Lambda$ができます。(ラムダの大文字です。)

固有値で表現することで$A$の特徴が見えやすくなったり、連続して積をとると左右に$P^{-1}P=I$が続くので$A^{n}$が計算しやすくなったりするというメリットがあります。なんか高校数学でやった気がする。

特に前者の目的がデータサイエンスでは重要。実際のデータ(特徴量とか特徴行列とか)は正方行列でないことがほとんどなので、より一般化した「特異値分解」という問題で扱うことになります。このまま読み進んでいけば、多分割とすんなり理解できる形で最後の方に出てきます。

で、上の式の$P$は何者かというと、直交行列です。

直交行列 orthogonal matrix

転置行列と逆行列が等しくなる正方行列のこと。下の関係を満たす。

つまり転置と元の行列の積が、どっちからかけても単位行列になる。

どんなやつかっていうと、たとえばこんなやつ。

NumPyで実際に確認してみる。

R = np.sqrt(1/2) * np.array([[1,1],[-1,1]])
R.T - np.linalg.inv(R) # 逆行列と同じか？
# array([[ 1.11022302e-16, -1.11022302e-16],
#        [ 1.11022302e-16,  1.11022302e-16]])
R@R.T # 転置と元をかけて、単位行列になるか？
# array([[ 1.00000000e+00, -4.26642159e-17],
#        [-4.26642159e-17,  1.00000000e+00]])

まあ、誤差があるので厳密には0になりませんでしたが、ほとんど0ということで、関係が成り立ってそうですね。

また、行列式が1か-1になるという性質もあります。

np.linalg.det(R) # 1.0

特異値 singular value と特異ベクトル singular vector

やっと本題が見えてきました。

特異値$\sigma$は、固有値$\lambda$を正方行列以外の$m \times n$行列に一般化したものです。

今まで$A$を使っていたので、紛れないようにこの$m \times n$行列を、行列$M$とします。以下の関係が成り立ちます。

$\vec{v}$と$\vec{u}$はそれぞれ特異ベクトルというもので、それぞれ$MM^{T}$と$M^{T}M$の固有ベクトルに相当します。($MM^{T}$は$m \times m$、$M^{T}M$は$n \times n$の正方行列なので、固有ベクトルが定義できます。)

ん、でも固有値って線型変換をスカラー倍にできるってことじゃなかったっけ...なんで$\vec{v}$と$\vec{u}$って別のベクトルになっとるんだ？ってのは一旦置いておきましょう。大事なのは任意の行列を次の対角行列に分解ができること。

特異値分解 singular value decomposition, SVD

正方行列ではない$m \times n$の行列$M$について、固有値分解のように、次の変形をすることを特異値分解という。

$U$は$m \times m$のユニタリ行列(unitary matrix)で、複素数の直交行列に相当する行列。普通は実数データしか考えないとして、ただの直交行列と思えばOK。同様に$V$は$n \times n$のユニタリ行列ですが、実数で考えて直交行列。

$S$は$M$の特異値$\sigma$を対角成分に持つ対角行列。ただし、対角行列である(正方行列の形で対角成分をもつ)には、行と列の形が一致する($m=n$である)必要があります。より一般化して考えると、$m>n$ならその対角行列の下に零行列で埋めて、$n>m$なら右側が零行列で埋めます。いずれにせよ、特異値(と零)が各列・各行に1つずつ並んだ行列になるので、行列$M$の特徴が現れた行列、と解釈できるようになります。

で、$U$の中身は上述の特異ベクトル$\vec{u}$を並べたもの、$V$の中身は特異ベクトル$\vec{v}$を並べたものだったりします。特異値分解したときに特異値の左右に来るので、それぞれ$\vec{u}$を左特異ベクトル、$\vec{v}$を右特異ベクトルと呼ぶようです。

NumPyで書くとこうなります。

M = np.array([[ 1,  2,  3,  4],
              [ 5,  6,  7,  8],
              [ 9, 10, 11, 12]]) # m:3, n:4の3行4列。

np.linalg.matrix_rank(M) # 2 : ランクがmより小さい必要あり。

u,s,vh = np.linalg.svd(M,full_matrices=False)
# SVDが特異値分解。
# u,s,vhにそれぞれ左特異ベクトルの行列、特異値、右特異ベクトルの行列を返す。
# full_matrices=Trueならuとvhのサイズは(m,m)、(n,n)
# full_matrices=Falseならuとvhのサイズは(m,k)、(k,n)、k=min(m,n)

# 一応サイズを確認して、ドット積が使えることを確認。
print(u.shape, s.shape, vh.shape) # (3, 3) (3,) (3, 4)

# 中身を確認
u # これが左特異ベクトルを並べた行列
# array([[-0.20673589, -0.88915331,  0.40824829],
#        [-0.51828874, -0.25438183, -0.81649658],
#        [-0.82984158,  0.38038964,  0.40824829]])

s # これが行列Mの特異値
# array([2.54368356e+01, 1.72261225e+00, 5.14037515e-16])

vh # これが右特異ベクトルを並べた行列(の転置)
# array([[-0.40361757, -0.46474413, -0.52587069, -0.58699725],
#        [ 0.73286619,  0.28984978, -0.15316664, -0.59618305],
#        [ 0.44527162, -0.83143156,  0.32704826,  0.05911168]])

# ドット積で元のMと同じになるか確認。sはdiagで対角行列Sにして計算する。
u@np.diag(s)@vh
# array([[ 1.,  2.,  3.,  4.],
#        [ 5.,  6.,  7.,  8.],
#        [ 9., 10., 11., 12.]]) : Mと同じになった!

というわけで、無事に任意の行列$M$について特異値と特異ベクトルを求めることができ、それぞれが何者なのか大体理解頂けたと思います。

参考 : NumPy | numpy.linalg.svd

エルミート行列 hermitian transpose、随伴行列 adjoint matrix

ところで、なぜNumPyの特異値分解の3つめの返り値は転置(Transpose)のTをつかったvtじゃなくてvhなのかと気になりませんでしたでしょうか？ならない

ドキュメントでも転置が$A^{T}$ではなく$A^{H}$と書かれているのは、特異値分解というものが複素数にも対応しているためで、エルミート転置(Hermitian transpose)の記号が使われているためです。これは、転置したあと、複素共軛、つまり全成分の虚部の符号を全て逆転したもので、随伴行列(adjoint matrix)とも言います。

とはいえ複素数のデータを扱わない限りお世話にならないので、ここではエルミート行列の詳細は省きますが、要は、実数のデータを扱う場合は、エルミート行列=転置した行列と読み替えてもらって問題ないなので、hでもtでも良さそうということです。

こういう言葉知ってるだけでも、ドキュメントがうんと読みやすくなりますね。

むすび

と、いうわけでデータサイエンスに必要な線型代数の知識が一瞬(?)で復習できました。思ったよりボリュームが多くなってしまいましたが、理系の方の復習には良い感じになってると思います。

特に今回は統計やデータサイエンスに頻出のテクニックである主成分分析に必要な「特異値」「特異ベクトル」「特異値分解」までの理解を念頭に、要点を押さえやすい構成を意識して書いてみました。これでデータサイエンス入門記事の主成分分析編が書ける、はず...。

何かまた必要になったら、随時書き足していこうと思います。

それでは〜