ノイズ掃き出し法 with Python

はじめに

今回は、高次元主成分分析(高次元PCA)の1つである、ノイズ掃き出し法(noise-reduction 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)}$ とする。

標本共分散行列$S_n$の固有値分解は

$$\begin{align*} S_n&=\frac{1}{n-1}(X-\bar X)(X-\bar X)^T\\ &=\sum_{s=1}^d\hat\lambda_{(s)}\hat h_{(s)}\hat h_{(s)}^T \end{align*}$$

双対標本共分散行列$S_{D,n}$の固有値分解は

$$\begin{align*} S_{D,n}&=\frac{1}{n-1}(X-\bar X)^T(X-\bar X)\\ &= \sum_{s=1}^n\hat \lambda_{(s)}\hat u_{(s)}\hat u_{(s)}^T \end{align*}$$

$S_n$の最初のn個の固有値$\lambda_{(1)} \ge \cdots \ge \lambda_{(n)}(\ge0)$は$S_{D,n}$の固有にもなっている。

一般的なPCAは標本共分散行列をもちいて$\Sigma$を推定する。ここでは次元$d>n$なので計算量を減らすために双対標本共分散行列を用いている。やっていることはPCAと同じ。

$X-\bar X$の特異値分解は$S_n$と$S_{D,n}$の固有値、固有ベクトルを用いて

$$\frac{X-\bar X}{\sqrt{n-1}}=\sum_{s=1}^{\min\{d,n-1\}} \hat\lambda_{(s)}^{1/2}\hat h_{(s)} \hat u_{(s)}^T$$

これより

$$\hat h_{(s)}=\frac{X-\bar X}{\sqrt{(n-1)\hat\lambda_{(s)}}}\hat u_{(s)}\hspace{5pt}(s=1,…,\min\{d,n-1\} )$$

ノイズ掃き出し法

ノイズ掃き出し法による固有値は

$$\tilde\lambda_{(s)}=\hat\lambda_{(s)}-\frac{tr(S_{D,n})-\sum_{t=1}^s\hat\lambda_{(t)}}{n-1-s} \hspace{5pt} (s=1,...,\min\{d,n-2\})$$

ノイズ掃き出し法による固有ベクトルは

$$\tilde h_{(s)}=\frac{(X-\bar X)}{\sqrt{(n-1)\tilde \lambda_{(s)}}}\hat u_{(s)}=\sqrt{\frac{\hat \lambda_{(s)}}{\tilde\lambda_{(s)}}}\hat h_{(s)} \hspace{5pt}(s=1,...,\min\{d,n-2\})$$

ノイズ掃き出し法による主成分スコア

$$\tilde s_{(k)}=\hat u_{(k)}\sqrt{(n-1)\tilde \lambda_{(k)}}\hspace{5pt} (k=1,...,\min\{d,n-2\})$$

Pythonコード

import numpy as np

#noise reduction methodology
def NR(X):
    d, n = X.shape

    # 平均
    one_n = np.ones((n, 1))
    X_bar = X @ (one_n @ one_n.T) / n
    # 標本共分散行列
    S = (1/(n-1)) * (X - X_bar) @ (X - X_bar).T
    # 双対標本共分散行列
    S_D = (1/(n-1)) * (X - X_bar).T @ (X - X_bar)

    # S_Dの固有値分解
    val, vec = np.linalg.eigh(S_D)
    index = np.argsort(val)[::-1]
    lam_hat = val[index]
    u_hat = vec[:, index]

    # h_hatを求める
    h_hat = [(X - X_bar) @ u_hat[:, i] / np.sqrt((n-1) * lam_hat[i]) for i in range(min(d,n-2))]

    #S_Dのトレース
    tr_S_D = np.trace(S_D)

    # 固有値
    lam_tilde = [lam_hat[i] - (tr_S_D - np.sum(lam_hat[:i+1])) / (n - 2 - i) for i in range(min(d,n-2))]
    # 固有ベクトル
    h_tilde = [np.sqrt(lam_hat[i] / lam_tilde[i]) * h_hat[i] for i in range(min(d,n-2))]
    # 主成分スコア
    s_tilde = [u_hat[:, i] * np.sqrt((n-1) * lam_tilde[i]) for i in range(min(d,n-2))]

    return lam_tilde, h_tilde, s_tilde

まとめ

高次元データは扱うのが大変ですが、この方法によってかなり扱えるようになりました。
詳しく知りたい方は、参考文献にある高次元の統計学を見てみてください。

参考文献

• 青嶋誠、矢田和善 (2019)「高次元の統計学」
  https://www.kyoritsu-pub.co.jp/book/b10003167.html