Zenn
🟠

クロスデータ行列法 with Python

2023/11/02に公開
Python
math
統計
統計学
データサイエンス
tech

はじめに

今回は、高次元主成分分析(高次元PCA)の1つである、クロスデータ行列法(cross-data-matrix methodology)を紹介します。

n>dのPCAはこちらの記事をチェック。
https://zenn.dev/totopironote/articles/aa17833ef00e5f

高次元データについて知りたい方はこちらをチェック
https://zenn.dev/totopironote/articles/f13fe705f52d42

高次元データと主成分分析について知りたい方はこちらをチェック
https://zenn.dev/totopironote/articles/16807d3a809ad7

基本

d>nの状況を考える。

平均\mu , 共分散行列\Sigma をもつ母集団からn個のd次データベクトルx_j , j=1,…,nを無作為に抽出する。
d \times n のデータ行列をX = (x_1,...,x_n)と定義する。x_j=(x_{j(1)},...,x_{j(d)})^T

\Sigma の固有値分解を次のようにする。

\Sigma=\sum_{s=1}^d\lambda_{(s)}h_{(s)}h_{(s)}^T

ここで、 \lambda_{(1)} \ge \cdots \ge \lambda_{(d)} とする。

クロスデータ行列法

Step1.

データ行列XをX_l=(x_{1,l},\dots,x_{n(l),l}) ,l=1,2 に分割する。

n_{(1)}= \lceil n /2 \rceil , n_{(2)}=n-n_{(l)}

\bar X_l = (\bar x_l,\dots,\bar x_l)=\bar x_l1_{n(l)}^T ,\hspace{5pt} \bar x_l = n_{(l)}^{-1} \sum_{j=1}^{n_{(l)}} x_{j,l} \hspace{5pt} (l=1,2)

Step2.

クロスデータ行列を計算する。

\begin{align*} S_{D(1),n}&=\frac{(X_1-\bar X_1)^T(X_2-\bar X_2)}{\sqrt{(n_{(1)}-1)(n_{(2)}-1)}}\\ S_{D(2),n}&=S_{D(1),n}^T \end{align*}

Step3.

対称行列S_{D(1),n}S_{D(2),n} の固有値 \acute\lambda_{(1)}^2 \ge \cdots\ge \acute\lambda_{(n_{(2)}-1)}^2 \ge 0 と対応する固有ベクトル\acute u_{(s),1} を計算する。\lambda_{(s)}を\acute \lambda_{(s)} = \sqrt{\acute \lambda_{(s)}^2} で推定する。

Step4.

対称行列S_{D(2),n}S_{D(1),n} の固有ベクトル\acute u_{(s),2} を計算し、

\acute u_{(s),2} =Sign(\acute u_{(s),1}^T S_{D(1),n} \acute u_{(s),2})\acute u_{(s),2} で符号を調整する。sごとに調整していることに注意。

Step5.

\acute h_{(s)} = \frac{1}{\sqrt {\acute \lambda_{(s)}}}\left( \frac{(X_1 - \bar X_1)\acute u_{(s),1}}{\sqrt{n_{(1)}-1}}+ \frac{(X_2 - \bar X_2)\acute u_{(s),2}}{\sqrt{n_{(2)}-1}} \right)

を計算し、固有ベクトルh_{(s)}を\acute h_{(s)*}=\acute h_{(s)}/\|\acute h_{(s)}\|で推定する。

Step6.

\acute u_{(s),l}=(\acute u_{1(s),l}, \dots ,\acute u_{n_{(l)}(s),l})^T \hspace{5pt} (l=1,2) とする。

第k主成分スコアは

\acute s_{j(k),l} = \acute u_{j(k),l}\sqrt{(n_{(l)}-1)\acute \lambda_{(k)}} \hspace{5pt} (l=1,2)

そして各kで\acute s_{j(k)}=\acute s_{1(k),1},\dots \acute s_{n_{(1)}(k),1},\acute s_{1(k),2},\dots \acute s_{n_{(2)}(k),2}\hspace{5pt} (j=1,\dots ,n) とする。

いまクロスデータ行列法で固有値などを推定する方法を紹介した。この手順は高次元の統計学に載っている。Step3,4で対称行列にしてから固有値分解を行なっている。

ここで以下の特異値分解ができればよいのでは？という疑問が生まれた。

S_{D(1),n} = \sum_{s=1}^{n_{(2)}-1}\acute \lambda_{(s)} \acute u_{(s),1}\acute u_{(s),2}^T

実際 Pythonでは特異値分解を計算できるので、やってみるとちょっとズレる。

Pythonコード

高次元の統計学で用いられているやり方。固有値分解を使用。

#cross-data-matrix methodology 固有値分解

def CDM1(X):
    # Step1
    d, n = X.shape
    n1 = np.ceil(n/2).astype(int)
    n2 = n - n1
    X1 = X[:, :n1]
    X2 = X[:, n1:]

    # Step2
    X1_bar = np.mean(X1, axis=1).reshape(-1, 1)
    X2_bar = np.mean(X2, axis=1).reshape(-1, 1)
    S_D1 = (X1 - X1_bar).T @ (X2 - X2_bar) / np.sqrt((n1 - 1) * (n2 - 1))
    S_D2 = S_D1.T

    # Step3 S_D1@S_D2の固有値・固有ベクトル
    l1, u1 = np.linalg.eigh(S_D1@S_D2)
    index = np.argsort(l1)[::-1][:(n2-1)]
    acute_lam = np.sqrt(l1[index])
    acute_u1 =u1[:, index]

    # Step4 S_D2@S_D1の固有値・固有ベクトル
    l2, u2 = np.linalg.eigh(S_D2@S_D1)
    index2 = np.argsort(l2)[::-1][:(n2-1)]
    u2 = u2[:, index2]
    acute_u2 =np.array([np.sign(acute_u1[:,i].T @ S_D1 @ u2[:,i]) * u2[:,i] for i in range(n2-1)]).T

    # Step5　固有ベクトル
    acute_h1 = np.array([(X1 - X1_bar) @ acute_u1[:, i]/np.sqrt((n1 - 1) * lam) for i, lam in enumerate(acute_lam)])
    acute_h2 = np.array([(X2 - X2_bar) @ acute_u2[:, i]/np.sqrt((n2 - 1) * lam) for i, lam in enumerate(acute_lam)])
    acute_h = np.array([(h1 + h2) /np.linalg.norm(h1+h2) for h1, h2 in zip(acute_h1, acute_h2)])

    # Step6 主成分スコア
    s1 = acute_u1 * np.sqrt((n1 - 1) * acute_lam)
    s2 = acute_u2 * np.sqrt((n2 - 1) * acute_lam)
    s = np.vstack((s1,s2))

    return acute_lam, acute_h, s

特異値分解を使用するやり方。

#cross-data-matrix methodology　特異値分解

def CDM2(X):
    # Step1
    d ,n = X.shape
    n1 = np.ceil(n/2).astype(int)
    n2 = n - n1
    X1 = X[:, :n1]
    X2 = X[:, n1:]

    # Step2
    X1_bar = np.mean(X1, axis=1).reshape(-1, 1)
    X2_bar = np.mean(X2, axis=1).reshape(-1, 1)
    S_D1 = (X1 - X1_bar).T @ (X2 - X2_bar) / np.sqrt((n1 - 1) * (n2 - 1))
    S_D2 = S_D1.T

    # Step3 , Step4　固有値
    u1, l, u2 = np.linalg.svd(S_D1, full_matrices=False)
    acute_lam = l[:(n2-1)]
    acute_u1 = u1[:,:(n2-1)]
    acute_u2 = u2[:(n2-1)].T
    

    # Step5　固有ベクトル
    acute_h1 = np.array([(X1 - X1_bar) @ acute_u1[:, i]/np.sqrt((n1 - 1) * lam) for i, lam in enumerate(acute_lam)])
    acute_h2 = np.array([(X2 - X2_bar) @ acute_u2[:, i]/np.sqrt((n2 - 1) * lam) for i, lam in enumerate(acute_lam)])
    acute_h = np.array([(h1 + h2) /np.linalg.norm(h1+h2) for h1, h2 in zip(acute_h1, acute_h2)])

    # Step6 主成分スコア
    s1 = acute_u1 * np.sqrt((n1 - 1) * acute_lam)
    s2 = acute_u2 * np.sqrt((n2 - 1) * acute_lam)
    s = np.vstack((s1,s2))

    return acute_lam, acute_h, s

CDM1は固有値分解を用いた。一方でCDM2は特異値分解を用いた。

CDM1とCDM2によって求められた固有値は一致した。しかし固有ベクトルの符号が一致しなかった。そのため主成分スコアも符号が一致しなかった。とりあえず、高次元の統計学にはCDM1が使われているのでCDM1をつかのが無難である。

まとめ

この方法はシンプルにすごいなと思いました。普通ならブートストラップみたいにデータ数を増やす方向を考えると思います。しかしデータを半分にして掛け合わせることで、小さいところはさらに小さく、大きいところはさらに大きくなるようにしています。すごい発想力。

参考文献

Discussion