多次元尺度構成法イントロダクション

ちょっと多次元尺度構成法について理解する必要があったので簡単なまとめです。

多次元尺度構成法とは、多数の多次元点間の距離データのみが与えられたときに、その距離を再現するような空間(座標系)を逆算する手法です。(多次元)座標値が与えられれば距離は自由に計算できますが、距離だけが与えられた場合に座標値を計算するのは直感的には中々できません。

今、n個の点がありそれぞれの点間の $n^2-n$ 個の距離 $d_{ij} \quad i,j\in \{1,\cdots,n\}$ がデータとして与えられたとします( $d_{ii}$ はゼロです)。空間の次元をqとし、未知のn個の座標ベクトルを $\bf{x}_i \in \mathcal{R}^q, \quad i \in {1,\cdots,n}$ 、それらを縦に並べた行列を

$\bf{X} = \left( \begin{array} \bf{x}_1^t \\ \vdots \\ \bf{x}_n^t \end{array}\right)$

とします(tは転地記号)。また、点間の距離はユークリッド距離

$d_{ij}^2 = (\bf{x_i} - \bf{x_j})^2 = \sum_k \left(x_{ik}x_{ik} + x_{jk}x_{jk} - 2*x_{ik}x_{jk}\right)$

とします。
ここで

$\bf{B} = \bf{X}\bf{X}^t$

という行列を考えます。もし何らかの形で距離データからBを求めることができれば、Bを対角化し固有値の平方根をとることで、 $\bf{X}$ を求めることができます。 $\bf{B}$ の成分は

$b_{ij} = \sum_k x_{ik}x^t_{kj} = \sum_kx_{ik}x_{jk}$ 　　　　　　　　(a)

であり、距離行列の成分と比べて次のような関係にあることがわかります。

$d_{ij}^2 = b_{ii} + b_{jj} - 2b_{ij}$ 　　　　　　　　(b)

(少々ややこしいですが、ここではi,jについては和を取りません。左辺にも現れるもので、ダミー変数ではないとわかります。)
この関係式によって、距離行列(手元にあるデータ)と $\bf{B}$ を結びつけることができます。

ところが、実はこの関係式だけでは距離行列から $\bf{B}$ を(すなわち座標値を)決めることができません。考えてみれば当然で、座標軸は原点をずらしたり回転したりしても距離は変わりません。逆に、距離だけがわかっている場合は原点や回転の自由度があるため座標軸を一意に決めることができないのです。
回転については直交変換の自由度が残ることになりますが、原点の位置を決めればまずは一つの座標系を書き下すことができるようになります。普通、原点は計算したデータ座標値の平均がゼロとなるようにとられます。

$\sum_{i=1}^n x_{ik} = 0, \quad k\in \{1,\cdots,q\}$

この関係式を(a)式に適用すると、 $b_{ij}$ もまたiまたはjの一方について和をとるとゼロとなることがわかります。これにより、(b)についてiで和をとると

$\sum_i d_{ij}^2 = \text{Tr}(\bf{B}) + nb_{jj}$

となります。ここで $\text{Tr}(\bf{B})$ は $\bf{B}$ の対角成分の和(トレース)で、nは $\bf{B}$ の次元です。同様に、jについての和、i,j両方についての和をとると

$\sum_j d_{ij}^2 = nb_{ii} + \text{Tr}(\bf{B})$

$\sum_{i,j} d_{ij}^2 = 2n\text{Tr}(\bf{B})$

という関係式も得られ、これら3つを $b_{ii},b_{jj},\text{Tr}(\bf{B})$ について解くことで、(b)式より $b_{ij}$ を $d_{ij}^2$ を用いて書き下すことができます。

$b_{ij} = -\frac{1}{2}\big( d_{ij}^2 - d_{i\cdot}^2 - d_{\cdot j}^2 + d_{\cdot \cdot}^2 \big)$
$d_{i\cdot}^2 = \frac{1}{n}\sum_j d_{ij}^2$
$d_{\cdot j}^2 = \frac{1}{n}\sum_i d_{ij}^2$
$d_{\cdot \cdot}^2 = \frac{1}{n^2}\sum_{ij} d_{ij}^2$
これで $\bf{B}$ がわかったので、ここから $\bf{B}=\bf{X}\bf{X}^t$ を解いて $\bf{X}$ を求めます。そのためには $\bf{X}$ の固有値・固有ベクトルを求めればよいというのが線形代数の常套手段です。実際

$\bf{B}\bf{x}^k = \lambda_k \bf{x}^k$

となる固有値 $\lambda_k$ 、固有ベクトル $\bf{x}^k$ が求まったとすると、

$\bf{B} = \sum_k\lambda_k\bf{x}^k\bf{x}^{kt}$

と書くことができ、これは、固有ベクトルを並べた行列

$\bf{V} = (\bf{x}_1, \cdots , \bf{x}_n)$

と、固有値を対角成分に並べた行列 $\bf{\Lambda} = \text{diag} \{\lambda_1, \cdots , \lambda_n\}$ を用いて、

$\bf{B} = \bf{V}\bf{\Lambda}\bf{V}^t = \bf{V}\bf{\Lambda}^{1/2}\left(\bf{V} \bf{\Lambda}^{1/2}\right)^t$

と書けます。ただし $\bf{\Lambda}^{1/2} = \text{diag} \{\sqrt{\lambda_1}, \cdots , \sqrt{\lambda_n}\}$ です。
従って、求める $\bf{X}$ は

$\bf{X} = \bf{V}\bf{\Lambda}^{1/2}$

であることがわかります。
ここで注意が必要です。固有値の平方根を要素とする行列を作りましたが、固有値は一般に負の数や複素数に成り得るので、この議論が成り立つためには固有値は0または正でなければなりません。これは行列 $\bf{B}$ が半正定値でなければならないということを意味します。実際には必ずしも半正定値ではなく、また欲しいのはなるべくズレの少ない低次元表現ですので、固有値を降順にソートし、大きいものからp個持ってきた行列で $\bf{\Lambda}^{1/2}$ を近似します。これによってp次元で表現された座標値 $\bf{X}$ が求まります。

長々と解説しましたが、Rでこれを実行するのは簡単です。statsパッケージにcmdscaleという関数があり、デフォルトでは2次元の多次元尺度構成法を行ってくれます。

> library(stats)
> iris.dist <- dist(iris[,-5])        #距離行列を作成
> iris.cmd <- cmdscale(iris.dist,k=2) #多次元尺度構成法
> plot(iris.cmd,col=iris$Species)     #結果を色つきでプロット

以下のようなグラフが得られます。

irisは元々4次元のデータですが、多次元尺度構成法によって得られた2次元表現でも、そこそこきちんと種別を反映してくれていることがわかると思います。irisは元データが座標値に対応するため、低次元表現を求めるという意味では主成分分析と同等ですが、データによってはそもそも距離データしか得られないということもあり、そういう場合には得にデータの可視化に役立ちます。

今回紹介した多次元尺度構成法は古典的多次元尺度構成法と呼ばれるもので、実際には距離がユークリッドでない場合や非線形な場合などを主眼に置いた様々な改良手法が存在します。いずれその辺りも書くかもしれません。