概要
データ分析について初めて触れる機会がありましたので、備忘録として計量的多次元尺度法について調べた内容を整理してみました。データ分析、pythonおよび数学については素人なので内容に誤りなどあるかもしれませんがご容赦ください。
多次元尺度法とは
与えられたデータから何らかの非類似度（距離など）を求め、2次元などの低次元での座標値を算出する手法です。これにより、多次元データであっても、2次元散布図などから相対的な位置関係を把握することができます。
例えば、「東京／大阪などの都市間の距離や移動時間」という非類似度のみを入力データとして、「相対的な位置関係」を得ることができます。都市間の距離は緯度経度などの低次元データから計算で求まるかと思いますが、これが任意の多次元であっても低次元に落とし込むことができ、データ間の大枠での関係性を見ることができます。
距離の公理
非類似度の尺度として使用する距離にはいくつか種類がありますが、以下の距離の公理を満たす必要があります。
| 名前 | 内容 | 数式 |
|---|---|---|
| 非負性 | 距離は正となる | $d_{ij} \geq 0$ |
| 非退化性 | 始点と終点が一致する場合は0になる | $d_{ii} = 0$ |
| 対称性 | 始点と終点を入れ替えても一致する | $d_{ij} = d_{ji}$ |
| 三角不等式が成立 | 任意の中点を経由せず結ぶ直線が最短となる | $d_{ij} \leq d_{ik} + d_{kj}$ |
計量的多次元尺度法の流れ
実行してる内容としてはざっくり以下となります。
（１）非類似度（距離）の2乗を要素とする行列を作成する
（２）ヤング・ハウスホルダー変換を行う
（３）固有値、固有ベクトルを算出する
（４）スペクトル分解を行い、固有値を近似する
（５）相対座標を取得
理論的なところ
とりあえず理論的なところを多少は理解しておきたいので、何年かぶりに数式に向き合ってみます。まず、前提として以下のような状態を考えてみます。
| 項目 | 表記 |
|---|---|
| 入力データの数 | $n$ |
| 各データの持つ要素（次元） | $r$（ただし、$r \lt n$と想定） |
| $i$番目の入力データの座標 | $(x_{i1}, x_{i2}, ... , x_{ir})$ |
| $i$番目の入力データの位置ベクトル | $X_i$ |
距離2乗行列の表現
$i$番目と$j$番目のデータ点間の距離（ユークリッド距離）の2乗は以下のようにして算出できます。
d_{ij}^2 = \sum_{k=1}^{r}(x_{ik} - x_{jk})^2
これを要素とする距離2乗行列を$D^2$とすると、以下のように$n×n$の正方行列となります。ここで、距離の公理の非退化性と対称性より、行列の対角成分が0である対象行列となっています。
\textbf{D}^2 =
\begin{pmatrix}
d_{11}^2 & d_{12}^2 & \cdots & d_{1n}^2 \\
d_{21}^2 & d_{22}^2 & \cdots & d_{2n}^2 \\
\vdots & \vdots & \ddots &\vdots \\
d_{n1}^2 & d_{n2}^2 & \cdots & d_{nn}^2
\end{pmatrix}
=
\begin{pmatrix}
0 & d_{12}^2 & \cdots & d_{1n}^2 \\
d_{21}^2 & 0 & \cdots & d_{2n}^2 \\
\vdots & \vdots & \ddots &\vdots \\
d_{n1}^2 & d_{n2}^2 & \cdots & 0
\end{pmatrix}
この行列$\textbf{D}^2$を、入力データの座標情報$(x_{i1}, x_{i2}, ... , x_{ir})$（ただし、$1 \leq i \leq r$）を各成分とする座標行列$\textbf{X}$を用いて表現すると、以下のように表すことができます。多次元尺度法では、この式①を使用していきます。
\textbf{D}^2 = diag(\textbf{XX}^T)\textbf{11}^T - 2\textbf{XX}^T + \textbf{11}^Tdiag(\textbf{XX}^T) \hspace{1cm} \cdots \:①
ただし、ここでは以下としています。
$\textbf{X}$は、行が1点のデータに対応する$n×r$行列。
\textbf{X} =
\begin{pmatrix}
x_{11} & \cdots & x_{1r} \\
\vdots & \ddots & \vdots \\
x_{n1} & \cdots & x_{nr} \\
\end{pmatrix}
$\textbf{11}^T$は、全成分が1の$n×n$正方行列。$\textbf{1}^T = \left(1, \cdots , 1 \right)$ に対応。
\textbf{11}^T =
\begin{pmatrix}
1 & \cdots & 1 \\
\vdots & \ddots & \vdots \\
1 & \cdots & 1 \\
\end{pmatrix}
$diag()$は対称行列の対角成分のみを抜き出したもの。ここでは、$\textbf{XX}^T$の対角成分のみを抜き出した$n×n$正方対角行列。
diag(\textbf{XX}^T) =
\begin{pmatrix}
\sum_{i=1}^r(x_{1i}^2) & 0 & \cdots & 0 \\
0 & \sum_{i=1}^r(x_{2i}^2) & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \sum_{i=1}^r(x_{ni}^2) \\
\end{pmatrix}
ヤング・ハウスホルダー変換
この行列$\textbf{D}^2$と$\textbf{X}$の関係式①に対して、ヤング・ハウスホルダー変換を行います。ヤング・ハウスホルダー変換では、式①の左右から中心化行列を作用させます。この操作により余分な項が落ち、中心化された座標行列のみで結果を表現することができます。
中心化行列は、$n$個の入力データの中心が原点に移るような操作に対応します。
データの平均座標を $\left( m_1, \ m_2, \ \cdots , \ m_r \right)$ とすると、以下のように表せます。
m_j = \frac{\sum_{i=1}^nx_{ij}}{n}
中心化された座標行列を$\textbf{X}^* $、その要素を$x_{ij}^* $とすると、中心化は$x_{ij}^* = x_{ij} - m_j $ に対応します。これを行列で書くと、
\begin{align}
\textbf{X}^* &=
\begin{pmatrix}
x_{11}^* & \cdots & x_{1r}^* \\
\vdots & \ddots & \vdots \\
x_{n1}^* & \cdots & x_{nr}^*
\end{pmatrix} \\
&=\begin{pmatrix}
x_{11} - m_1 & \cdots & x_{1r} - m_r \\
\vdots & \ddots & \vdots \\
x_{n1} - m_1 & \cdots & x_{nr} - m_r
\end{pmatrix} \\
&= \textbf{X}
- \begin{pmatrix}
m_1 & \cdots & m_r \\
\vdots & \ddots & \vdots \\
m_1 & \cdots & m_r
\end{pmatrix} \\
&= \textbf{X}
- \frac{1}{n}\begin{pmatrix}
\sum_{i=1}^nx_{i1} & \cdots & \sum_{i=1}^nx_{ir} \\
\vdots & \ddots & \vdots \\
\sum_{i=1}^nx_{i1} & \cdots & \sum_{i=1}^nx_{ir}
\end{pmatrix} \\
&= \textbf{X}
- \frac{1}{n}
\begin{pmatrix}
1 & \cdots & 1 \\
\vdots & \ddots & \vdots \\
1 & \cdots & 1
\end{pmatrix}
\begin{pmatrix}
x_{11} & \cdots & x_{1r} \\
\vdots & \ddots & \vdots \\
x_{n1} & \cdots & x_{nr}
\end{pmatrix} \\
&=\left(\textbf{I}-\frac{\textbf{11}^T}{n} \right) \textbf{X}
\\
\end{align}
ただし、$\textbf{I}$は$n×n$の単位行列です。
ということで、中心化行列を$\textbf{Q}$とすると、
\textbf{Q} = \textbf{I}-\frac{\textbf{11}^T}{n} \hspace{1cm} \cdots \: ②
となることがわかりました。
この中心化行列$\textbf{Q}$を式①の両側から作用させますが、$\textbf{Q1}=\textbf{1}^T\textbf{Q}=0$となるために右辺は第2項目のみが残ります。
\begin{align}
\textbf{QD}^2\textbf{Q} &=
\textbf{Q}\left(
diag(\textbf{XX}^T)\textbf{11}*^T - 2\textbf{XX}^T + \textbf{11}^Tdiag(\textbf{XX}^T)
\right)\textbf{Q} \\
&= -2\textbf{QXX}^T\textbf{Q} \\
&= -2\left(\textbf{QX}\right) \left(\textbf{QX}\right)^T \\
&= -2\textbf{X}^* \textbf{X}^{*T}
\end{align}
したがって、中心化された座標行列$\textbf{X}^* $について
\textbf{P} = \textbf{X}^* \textbf{X}^{*T} = -\frac{1}{2}\textbf{QD}^2\textbf{Q}
となります。行列$\textbf{P}$は式①と②から算出することができるため、具体的な値を求めることができます。
なお、以下よりPは対称行列となることがわかります。
\begin{align}
\textbf{P}^T &= \left(\textbf{X}^{*}\textbf{X}^{*T}\right)^T \\
&=\left(\textbf{X}^{*T} \right)^T \left(\textbf{X}^{*}\right)^T \\
&=\textbf{X}^{*}\textbf{X}^{*T} \\
&=\textbf{P}
\end{align}
スペクトル分解
ヤング・ハウスホルダー変換により得られた行列$\textbf{P}$は$n×n$の実対称行列となっているので、固有ベクトルからなる直行行列$\Gamma$と固有値のみを対角成分として持つ対角行列$\Lambda$を用いて以下のように表せます。
\begin{align}
\textbf{P} &=\textbf{X}^{*}\textbf{X}^{*T} \\
&=\Gamma\Lambda\Gamma^T \\
&=\Gamma\Lambda^{1/2}\Lambda^{1/2}\Gamma^T \\
&=\Gamma\Lambda^{1/2}\left(\Gamma\Lambda^{1/2}\right)^T
\end{align}
これより、中心化された座標行列について
\textbf{X}^{*} = \Gamma\Lambda^{1/2} \hspace{1cm} \cdots \: ③
となるので、固有値と固有ベクトルから求められることになります。なお、それぞれの行列については$i$番目の固有値を$\lambda_i$、それに対応する固有ベクトルを$\gamma_{i} = \left( \gamma_{1i} , \ \gamma_{2i} , \ \cdots , \ \gamma_{ni}\right)^T$とすると以下のようになっています。
\Lambda =
\begin{pmatrix}
\lambda_1 & & 0 \\
& \ddots & \\
0 & & \lambda_n
\end{pmatrix} \\
\Gamma =
\begin{pmatrix}
\gamma_{1} & \cdots & \gamma_{n}
\end{pmatrix}
=
\begin{pmatrix}
\gamma_{11} & \cdots & \gamma_{1n} \\
\vdots & \ddots & \vdots \\
\gamma_{n1} & \cdots & \gamma_{nn} \\
\end{pmatrix} \\
スペクトル分解では固有値について寄与率の低いものを近似的にゼロとみなすことで、出力となる座標情報の次元数を削減します。例えば、2次元の相対位置関係を得たい場合には$\lambda_1$と$\lambda_2$までを考慮し、$\lambda_3$以降の固有値はゼロとします。ここではスペクトル分解として、固有値を $m \left( \leq r\right)$ まで使用し、それ以降を0と近似するとします。この時、$\textbf{X}^{*}$は式③より以下となります。
\begin{align}
\textbf{X}^{*}
&= \Gamma\Lambda^{1/2} \\
&=
\begin{pmatrix}
\gamma_{1} & \cdots & \gamma_{n}
\end{pmatrix}
\begin{pmatrix}
\sqrt{\lambda_1} & & 0 \\
& \ddots & \\
0 & & \sqrt{\lambda_n}
\end{pmatrix} \\
&=
\begin{pmatrix}
\sqrt{\lambda_1}\gamma_{1} & \cdots & \sqrt{\lambda_m}\gamma_{m} & \sqrt{\lambda_{m+1}}\gamma_{m+1} & \cdots & \sqrt{\lambda_n}\gamma_{n}
\end{pmatrix} \\
&\sim
\begin{pmatrix}
\sqrt{\lambda_1}\gamma_{1} & \cdots & \sqrt{\lambda_m}\gamma_{m} & 0 & \cdots & 0
\end{pmatrix} \\
&=
\begin{pmatrix}
\sqrt{\lambda_1}\gamma_{11} & \cdots & \sqrt{\lambda_m}\gamma_{1m} & \\
\vdots & \ddots & \vdots & 0\\
\sqrt{\lambda_1}\gamma_{n1} & \cdots & \sqrt{\lambda_m}\gamma_{nm} &
\end{pmatrix} \hspace{1cm} \cdots \: ④
\end{align}
式④から、中心化された座標行列$\textbf{X}^{*}$の次元数$r$が、スペクトル分解により$m$次元まで削減されたことになります。これが多次元尺度法における出力値となります。
pythonで書いてみる
データ分析と言えばpythonだろう、ということで初めて触りながら書いてみました。sklearnを使えば簡単に実装できると思いますが、ここではあえて使わず、流れをつかむために以下に沿って進めてみました。
（１）非類似度（距離）の2乗を要素とする行列を作成する
（２）ヤング・ハウスホルダー変換を行う
（３）固有値、固有ベクトルを算出する
（４）スペクトル分解を行い、固有値を近似する
（５）相対座標を取得
（６）プロット
###（１）非類似度（距離）の2乗を要素とする行列を作成する
使用した入力データは放送大学の講義：データの分析と知識発見のこちらのものです。中身としては山手線の駅間の所要時間で構成された行列となっているため、今回はこれを距離とみなして進めます。
# -*- coding=utf_8_sig -*-
# データの読み込み
import pandas as pd
url="http://www.is.ouj.ac.jp/lec/16data/slide/Data/yamatesjis.csv"
df = pd.read_csv(filepath_or_buffer=url, encoding="shift_jis", sep=",", index_col=0)
# 距離2乗の行列を作成
d2 = df*df
###（２）ヤング・ハウスホルダー変換を行う
ここでは中心化行列を作成し、距離の2乗行列（d2）に作用させます。
import numpy as np
# 単位行列
e = np.eye(len(d2))
# 中心化行列
q = e - 1/len(d2)
# ヤングハウスホルダー変換
p = np.dot(np.dot( q, d2 ), q) * (-1/2)
###（３）固有値、固有ベクトルを算出する
import numpy.linalg as la
# 固有値の算出、固有ベクトルの算出
eig_val, eig_vec = la.eig(p)
###（４）スペクトル分解を行い、固有値を近似する
スペクトル分解で採用する固有値の数（出力で得られる座標の次元数）は変数dimで設定した2としています。また、出力となる位置関係の有意性の確認として寄与率の算出も行っています。
# 寄与率の算出
con_rate = []
tot_eig_val = np.sum(np.abs(eig_val))
for i in range(len(eig_val)):
con_rate.append(np.abs(eig_val[i])/tot_eig_val)
# スペクトル分解
dim=2
for i in range(dim):
print('Eigenvalue#' + str(i+1) + ' : ' + str(eig_val[i].real))
print('Contribution Rate = ' + str(con_rate[i]) + '\n----')
# 座標取得のため、固有値のルートを算出して固有ベクトルと積を取る
mod_eig_vec = np.sqrt(eig_val[i]) * eig_vec[:,i:i+1]
# 座標の取得
if i == 0:
crd = (mod_eig_vec)
else:
crd = np.append(crd, mod_eig_vec, axis=1)
###（５）相対座標を取得
得られた座標情報をファイルresult_mds.csvに書き出しておきます。
# 出力ファイル
outfile = open("result_mds.csv", encoding="utf_8_sig",mode="w")
# 結果の書き出し
i = 0
while i < len(df.iloc[1:1,1:].columns):
outfile.write(df.iloc[1:1,1:].columns[i])
j = 0
while j < dim:
outfile.write(',' + str(crd[i,j].real))
j+=1
i+=1
outfile.write('\n')
（６）プロット
得られた相対座標をmatplotlibで散布図として描画してみます。
import matplotlib.pyplot as plt
import japanize_matplotlib
import numpy as np
# ファイルの読み込み
dat = np.loadtxt("result_mds.csv", encoding="utf_8_sig", delimiter=",",dtype="unicode")
rows = dat.shape[0]
# グラフとフォントサイズの設定
plt.rcParams["font.size"] = 16
plt.figure(figsize=(8.0, 6.0))
# データの読み込み
label = dat[:,0:1]
x = dat[:,1:2].astype(np.float64)
y = dat[:,2:3].astype(np.float64)
# 座標軸にラベルを付ける
plt.xlabel("X")
plt.ylabel("Y")
# 散布図の作成
plt.scatter(x, y, c="b")
# ラベルの設定
for i in range(rows):
plt.annotate(label[i][0], xy=(x[i][0],y[i][0]))
# 散布図の描画
plt.show()
結果としては以下のようなものが得られ、距離情報から相対的な位置関係が得られたことになります。
最後に
基本的にはWebの情報中心で進めましたが、何となく多次元尺度法の流れや概要について理解ができたような気がします。初めてのPythonも楽しかったので、これをきっかけいろいろ手を出してみたいなーと思います。