🟠

POETとは？

2024/01/10に公開

Python

はじめに

POET(Principal Orthogonal complEment Thresholding)について解説します。POETは”詩人”ではないので注意です。笑

簡単に説明すると、高次元共分散行列の推定において疎な共分散行列分離させてその疎な共分散行列に対して閾値処理をすることで精度の高い推定をしています。

POET

共分散行列 $\Sigma$ の推定を考える。

まずデータ行列を $X$ ( $d\times n$ 行列)とする。dが次元数、nがデータ数である。このとき標本共分散行列は $S = (X-\bar X )(X - \bar X)^T/(n-1)$ となる。 $S$ の固有値を $\hat\lambda_1 \geq\dots \geq\hat\lambda_d$ とする。固有値に対応する固有ベクトルを $\hat h_1,\dots,\hat h_d$ とする。このとき $\Sigma$ の推定量を次で与える。

\hat \Sigma_K = \sum_{i=1}^K\hat \lambda_i\hat h_i\hat h_i^T + \hat R_K^\tau

とかける。ここで $\hat R_k = \sum_{i=k+1}^d\hat \lambda_i\hat h_i\hat h_i^T = (\hat r_{ij})_{d\times d}$ である。

$\Sigma$ の最初のk個の固有値が発散すると考える。kは未知。 $\hat R_K$ を次で定義する。

\hat{R}_K^\tau = (\hat{r}_{ij}^\tau)_{d \times d}, \quad \hat{r}_{ij}^\tau = \begin{cases} \hat{r}_{ii}, & \text{if } i = j; \\s_{ij}(\hat{r}_{ij})I(|\hat{r}_{ij}| \geq \tau_{ij}), & \text{if } i \neq j.\end{cases}

ここで $s_{ij}$ は一般化された収縮関数で、 $τ_{ij}>0$ は要素ごとに依存するしきい値である。

閾値処理はいくつもありますが、適応型閾値処理を考える。

\tau_{ij} = C \sqrt{\hat{\theta}_{ij}\omega}, \quad \omega = \frac{1}{\sqrt{d}} + \sqrt{\frac{\log d}{n}}

ここで $C>0$ はユーザーが決める。Cは推定量が正定値性を持つように選ばれる。

C_{min} = \text{inf}\{C > 0 : \lambda_{\min}(\hat{\Sigma}_{u,\hat K}^\tau(M)) > 0, \forall M > C\}

Cが十分に大きいと対角化される。Cは多重交差検証で推定することもできる。 $\hat \Sigma_{u,\hat K}^\tau = \hat R_{\hat K}^\tau$ である。Kの推定量が $\hat K$ である。

\hat w_j = (\hat w_{1j}, \ldots, \hat w_{dj})^T = (I_d - \sum_{i=1}^{K} \hat h_i \hat h_i^T) x_j

とした時、

\hat{\theta}_{ij} = n^{-1}\sum_{s=1}^{n} \left (\hat{w}_{is}\hat{w}_{js} - n^{-1} \sum_{s=1}^{n} \hat{w}_{is}\hat{w}_{js} \right)^2

ちなみに以前高次元主成分分析で解説したノイズ掃き出し法やクロスデータ行列法を用いればより精度の高い推定量が得られる。

Pythonコード

def th(M, t):
    return np.multiply(M, (np.abs(M) >= t).astype(int))

def poet(X, m):
    p, n = X.shape # X : p × n
    C = 1 # この値は論文や状況に応じて調整すること

    # BB^Tの推定
    val, vec = np.linalg.eigh(np.cov(X))
    index = np.argsort(val)[::-1][:m]
    hat_lam = val[index]
    hat_h = vec[:, index]
    hat_BBT = hat_h @ np.diag(hat_lam) @ hat_h.T

    # w_j の計算
    hat_w = np.array([(np.eye(p) - hat_h @ hat_h.T) @ x_j for x_j in X.T]).T

    # S_epsilon の推定
    hat_S_eps = hat_w @ hat_w.T / (n-1)

    # t_ij の計算
    hat_theta = np.zeros((p, p))  # 初期化
    sum = np.dot(hat_w, hat_w.T) / n  # すべてのペアワイズの和を計算
    squared_sum = np.outer(np.sum(hat_w, axis=1), np.sum(hat_w, axis=1)) / n**2  # 外積を利用して二乗和を計算

    hat_theta = sum - squared_sum  # 差を計算
    np.fill_diagonal(hat_theta, 0)  # 対角成分を0に設定
    t_ij = C * np.sqrt(hat_theta) * (p**(-1/2) + np.sqrt(np.log(p)/n))

    # Sigma の推定
    hat_Sigma = hat_BBT + th(hat_S_eps, t_ij)

    return hat_Sigma

まとめ

POETは疎な共分散行列を推定するのに適した推定法である。閾値処理で出てくるCの推定が厄介だと感じた。

参考文献

Jianqing Fan, Yuan Liao and Martina Mincheva. (2013) Large Covariance Estimation by Thresholding Principal Orthogonal Complements. https://arxiv.org/pdf/1201.0175.pdf