統計学実践ワークブックの第26章の内容について勉強したメモとなる。
多次元尺度法
$n$個体間の距離が与えられているとき、それら$n$個体の位置関係を座標で表現する手法として多次元尺度法（Multi-dimensional scaling；MDS）がある。多次元尺度法には、大きく分けて計量MDSと非計量MDSがある。
計量MDS
$n$個体間の距離データ$d_{ij}\geq 0,i\neq j$ $(i,j=1,2,\cdots,n)$が与えられている。ただし、$d_{ij}=d_{ji}$とする。$d_{ij}$の2乗を並べた行列（対角要素はすべて0とおく）$D$を距離行列とよぶ。
D=
\begin{pmatrix}
0 & d_{12}^2 & d_{13}^2 & \cdots & d_{1n}^2\\
d_{21}^2 & 0 & d_{23}^2 & \cdots & d_{2n}^2\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
d_{n1}^2 & d_{n2}^2 & d_{n3}^2 & \cdots & 0\\
\end{pmatrix}
この距離行列$D$をなるべく少ない次元数$k$の座標で表すことを考える。各個体の座標を表すベクトルを$\boldsymbol x_{i}=(x_{i1},x_{i2},\cdots,x_{ik})^\top,i=1,2,\cdots,n$とおく。ただし、各座標の平均値$\sum_{i=1}^{n}\boldsymbol x_{i}=\boldsymbol 0$とおく。このとき、$\sum_{i=1}^{n}|\boldsymbol x_i|^2=a$とおくと、
\begin{multline}
\begin{split}
d_{ij}^2 &= |\boldsymbol x_i-\boldsymbol x_j|^2=|\boldsymbol x_i|^2+|\boldsymbol x_j|^2-2\boldsymbol x_i^\top \boldsymbol x_j\\
\sum_{i=1}^n d_{ij}^2 &= \sum_{i=1}^n|\boldsymbol x_i-\boldsymbol x_j|^2=\sum_{i=1}^n(|\boldsymbol x_i|^2+|\boldsymbol x_j|^2-2\boldsymbol x_i^\top \boldsymbol x_j)=n|\boldsymbol x_i|^2+a\\
\sum_{i=1}^n\sum_{j=1}^nd_{ij}^2 &=\sum_{i=1}^n\sum_{j=1}^n|\boldsymbol x_i-\boldsymbol x_j|^2= \sum_{i=1}^n\sum_{j=1}^n(|\boldsymbol x_i|^2+|\boldsymbol x_j|^2-2\boldsymbol x_i^\top \boldsymbol x_j)=2na
\end{split}
\end{multline}
となり、
\begin{multline}
\begin{split}
\boldsymbol x_i^\top\boldsymbol x_j &= -\frac{1}{2}\left(d_{ij}^2-\frac{1}{n}\sum_{i=1}^nd_{ij}^2-\frac{1}{n}\sum_{j=1}^n d_{ij}^2+\frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^nd_{ij}^2 \right)\\
&=-\frac{1}{2}\left(1-\frac{1}{n}\sum_{i=1}^n\right)d_{ij}^2\left(1-\frac{1}{n}\sum_{j=1}^n\right)
\end{split}
\end{multline}
がなりたつ。$\boldsymbol X=(\boldsymbol x_1,\boldsymbol x_2,\cdots,\boldsymbol x_n)^\top$とおくと個体間の座標ベクトルの内積を並べた行列$\boldsymbol X\boldsymbol X^\top$は
\boldsymbol X\boldsymbol X^\top = -\frac{1}{2}\left(I_n-\frac{1}{n}J_n\right)D\left(I_n-\frac{1}{n}J_n\right)
となる。ここで、$I_n$は$n\times n$単位行列、$J_n$はすべての要素が１の$n\times n$行列。$\left(I_n-\frac{1}{n}J_n\right)$を左および右からかけるということは行列$D$の各行各列の平均を0にする変換であるため、二重中心化とよばれる。
$B=-\frac{1}{2}\left(I_n-\frac{1}{n}J_n\right)D\left(I_n-\frac{1}{n}J_n\right)$とおき、特異値分解を考える。$B$が非負定符号行列のとき、$B$の固有値$\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_r>0(r\leq n-1)$を並べた対角行列$\Lambda$、ノルム1に基準化された固有ベクトル$\boldsymbol u_1,\boldsymbol u_2,\cdots,\boldsymbol u_r$を並べた行列$U$より、$B=U\Lambda U^\top$となる。次元数$k=r$とすることで、$X=U\Lambda^{\frac{1}{2}}$となる。
このような特異値分解による行列の低ランク近似を計量心理学の分野ではエッカート・ヤング分解という。
$B$が非負定符号行列でないとき、どんなに次元数$k$を増やしてもユークリッド距離の2乗に基づいて距離行列$D$に一致する座標を求めることはできない。
次元数が大きすぎると結果を解釈するのが大変になるため、次元数を絞り込んで大きい方の固有値およびその固有ベクトルを用いる。
非計量MDS
与えられた距離データ$s_{ij}\geq0,i\neq j$$(i,j=1,2,\cdots,n)$を非計量値とみなして、大小関係のみをできるだけ保って単調変換を行う。
正準相関分析
2つの変数群$x_1,x_2,\cdots,x_p$および$y_1,y_2,\cdots,y_q$がある。このとき、それぞれの線形結合$a_1x_1+a_2x_2+\cdots +a_px_p$および$b_1y_1+b_2y_2+\cdots+b_qy_q$が最も高い相関を持つように係数$a_1,a_2,\cdots,a_p$、$b_1,b_2,\cdots,b_q$を定め、両変数群の関係を解釈していく手法を正準相関分析という。
変数群の標本分散共分散行列を$S_{xx}$および$S_{yy}$、変数群間の標本共分散行列を$S_{xy}$とおくと、$a_1x_1+a_2x_2+\cdots +a_px_p=\boldsymbol a^\top\boldsymbol x$および$b_1y_1+b_2y_2+\cdots+b_qy_q=\boldsymbol b^\top\boldsymbol y$の分散は、$\boldsymbol a^\top S_{xx}\boldsymbol a$および$\boldsymbol b^\top S_{yy}\boldsymbol b$となり、共分散は$\boldsymbol a^\top S_{xy}\boldsymbol b$となる。相関係数$\frac{\boldsymbol a^\top S_{xy}\boldsymbol b}{\sqrt{\boldsymbol a^\top S_{xx}\boldsymbol a\boldsymbol b^\top S_{yy}\boldsymbol b}}$を最大化すればよい。ただし、$\boldsymbol a$、$\boldsymbol b$の定数倍に対して相関係数は不変性があるため、
\max \boldsymbol a^\top S_{xy}\boldsymbol b \quad\text{s.t.}\quad \boldsymbol a^\top S_{xx}\boldsymbol a =1,\boldsymbol b^\top S_{yy}\boldsymbol b =1
を解く。
ラグランジュの未定乗数法より、
f(\boldsymbol a,\boldsymbol b,\lambda_1,\lambda_2)=\boldsymbol a^\top S_{xy}\boldsymbol b -\frac{\lambda_1}{2}(\boldsymbol a^\top S_{xx}\boldsymbol a -1)-\frac{\lambda_2}{2}(\boldsymbol b^\top S_{yy}\boldsymbol b -1)
$\boldsymbol a$、$\boldsymbol b$それぞれで偏微分して$\boldsymbol 0$とし、
\begin{multline}
\begin{split}
\frac{\partial f}{\partial \boldsymbol a}&=S_{xy}\boldsymbol b -\lambda_1 S_{xx}\boldsymbol a = \boldsymbol 0\\
\frac{\partial f}{\partial \boldsymbol b}&=S_{xy}^\top\boldsymbol a -\lambda_2 S_{yy}\boldsymbol b =\boldsymbol 0
\end{split}
\end{multline}
これより、$\lambda_1=\lambda_2=\boldsymbol a^\top S_{xy}\boldsymbol b$となり、$\boldsymbol a^\top S_{xy}\boldsymbol b=\lambda$、$\boldsymbol c = S_{xx}^{\frac{1}{2}}\boldsymbol a$として、
S_{xx}^{-\frac{1}{2}}S_{xy}S_{yy}^{-1}S_{xy}^\top S_{xx}^{-\frac{1}{2}}\boldsymbol c = \lambda^2\boldsymbol c
となる。これより、$\boldsymbol a^\top S_{xy}\boldsymbol b$が最大になるのは、$\lambda^2$が$S_{xx}^{-\frac{1}{2}}S_{xy}S_{yy}^{-1}S_{xy}^\top S_{xx}^{-\frac{1}{2}}$の第1固有値$\eta_1$、$\boldsymbol c$がその固有ベクトルとなる。また、$\boldsymbol a=S_{xx}^{-\frac{1}{2}}\boldsymbol c$、$\boldsymbol b = \frac{1}{\sqrt{\eta_1}}S_{yy}^{-1}S_{xy}^\top\boldsymbol a$となる。
例題
問26.1
[1]
\begin{multline}
\begin{split}
B&=-\frac{1}{2}\left(I_4-\frac{1}{4}J_4\right)D\left(I_4-\frac{1}{4}J_4\right)\\
&=-\frac{1}{2}\left(
\begin{pmatrix}
1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1
\end{pmatrix}
-\frac{1}{4}
\begin{pmatrix}
1&1&1&1\\
1&1&1&1\\
1&1&1&1\\
1&1&1&1
\end{pmatrix}
\right)
\begin{pmatrix}
0&5&5&16\\
5&0&4&5\\
5&4&0&5\\
16&5&5&0
\end{pmatrix}
\left(
\begin{pmatrix}
1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1
\end{pmatrix}
-\frac{1}{4}
\begin{pmatrix}
1&1&1&1\\
1&1&1&1\\
1&1&1&1\\
1&1&1&1
\end{pmatrix}
\right)\\
&=
\begin{pmatrix}
4&0&0&-4\\
0&1&-1&0\\
0&-1&1&0\\
-4&0&0&4
\end{pmatrix}
\end{split}
\end{multline}
$B$の固有値は掃き出し法をもちいればよい。
\begin{multline}
\begin{split}
\det
\begin{pmatrix}
4-\lambda&0&0&-4\\
0&1-\lambda&-1&0\\
0&-1&1-\lambda&0\\
-4&0&0&4-\lambda
\end{pmatrix}
&=(1-\lambda)(4-\lambda)
\begin{pmatrix}
1&0&0&-\frac{4}{4-\lambda}\\
0&1&-\frac{1}{1-\lambda}&0\\
0&-1&1-\lambda&0\\
-4&0&0&4-\lambda
\end{pmatrix}\\
&=(1-\lambda)(4-\lambda)
\begin{pmatrix}
1&0&0&-\frac{4}{4-\lambda}\\
0&1&-\frac{1}{1-\lambda}&0\\
0&0&\frac{\lambda^2-2\lambda}{1-\lambda}&0\\
0&0&0&\frac{\lambda^2-8\lambda}{4-\lambda}
\end{pmatrix}\\
&=(1-\lambda)(4-\lambda)\frac{\lambda^2-2\lambda}{1-\lambda}\frac{\lambda^2-8\lambda}{4-\lambda}
\begin{pmatrix}
1&0&0&-\frac{4}{4-\lambda}\\
0&1&-\frac{1}{1-\lambda}&0\\
0&0&1&0\\
0&0&0&1
\end{pmatrix}\\
&=\lambda^2(\lambda-2)(\lambda-8)
\end{split}
\end{multline}
よって固有値は、8、2、0、0となる。
その固有ベクトルも掃き出し法をもちいて計算すると(略)、$\left(\frac{1}{\sqrt{2}},0,0,-\frac{1}{\sqrt{2}}\right)^\top$、$\left(0,\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}},0\right)^\top$、$\left(\frac{1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\right)^\top$、$\left(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\right)^\top$となる。
また、Pythonで計算すると、
import numpy
B=numpy.array([[4,0,0,-4],[0,1,-1,0],[0,-1,1,0],[-4,0,0,4]])
B_eig=numpy.linalg.eig(B)
print(f"固有値は{B_eig[0]} 固有ベクトルは{B_eig[1]}")
固有値は[8. 0. 2. 0.] 固有ベクトルは[[ 0.70710678 0.70710678 0. 0. ]
[ 0. 0. -0.70710678 0.70710678]
[ 0. 0. 0.70710678 0.70710678]
[-0.70710678 0.70710678 0. 0. ]]
となり、一致する。
[2]
次元数1より、固有値8、固有ベクトル$\left(\frac{1}{\sqrt{2}},0,0,-\frac{1}{\sqrt{2}}\right)^\top$を用いる。
X=\left(\frac{1}{\sqrt{2}},0,0,-\frac{1}{\sqrt{2}}\right)*8^{\frac{1}{2}}=\left(2,0,0,-2\right)
$d_{ij}^2=|x_i-x_j|^2$なので、
D=
\begin{pmatrix}
0&4&4&16\\
4&0&0&4\\
4&0&0&4\\
16&4&4&0
\end{pmatrix}
[3]
次元数2より、固有値8、2、固有ベクトル$\left(\frac{1}{\sqrt{2}},0,0,-\frac{1}{\sqrt{2}}\right)^\top$、$\left(0,\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}},0\right)^\top$を用いる。
X=
\begin{pmatrix}
\frac{1}{\sqrt{2}}&0\\
0&\frac{1}{\sqrt{2}}\\
0&-\frac{1}{\sqrt{2}}\\
-\frac{1}{\sqrt{2}}&0
\end{pmatrix}
\begin{pmatrix}
8^{\frac{1}{2}}&0\\
0&2^{\frac{1}{2}}
\end{pmatrix}
=
\begin{pmatrix}
2&0\\
0&1\\
0&-1\\
-2&0
\end{pmatrix}
よって
D=
\begin{pmatrix}
0&5&5&16\\
5&0&4&5\\
5&4&0&5\\
16&5&5&0
\end{pmatrix}
となり、次元数を2に減らしたうえで、Dを完全に再現できる。