Zenn
🟠

ノイズ掃き出し法 with Python

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

はじめに

今回は、高次元主成分分析(高次元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

まとめ

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

参考文献

Discussion