PythonでPCA実装

この記事の目的

主成分分析（PCA）をPythonで実装してPCAがなにをやっているかを理解します。PCAの説明は他にも多数の解説記事があるので実装と関係ない部分の説明には踏み込みません。

主成分分析の理論

主成分分析では、データベクトル$\bf{x}$の分散が最大になる方向への線形変換を求める。

\bf{x}_i=\it{(x_{i1},...,x_{id})}^T\  \rm{(i=1,...,N)}

$N$個のデータからなるデータ行列を$\mathbf{X}=(\bf{x}_1,...,\bf{x}_N)^T$から平均ベクトル$\bf{\overline{x}}=\it{(\overline{x}_1,...,\overline{x}_d)}$を引いたデータ行列を$\bf{\overline{X}}$とする。

\bf{\overline{X}}=(x_1-\overline{x},...,x_N-\overline{x})^T \tag{1}

このデータ行列$\bf{\overline{X}}$の共分散行列$Var(\bf{\overline{X}})$の固有値$\lambda$が変換先成分の分散に対応し、固有ベクトル$\bf{a}_j$が変換ベクトルに対応する。

Var(\bf{\overline{X}})\bf{a}_j = \lambda\bf{a}_j \tag{2}

理論はこれだけなので、以下Pythonで実装する。

実装

共分散行列計算

  # 共分散行列計算
  ## N :データ数
  ## d :Xの特徴量次元数
  ## X.shape = (N, d)
  ## Sigma :共分散行列
  ## Sigma.shape = (d, d)
  N = X.shape[0]
  X_ = X - X.mean(axis=0).reshape(1,-1)
  Sigma = 1/N * X_.T @ X_

固有値、固有ベクトル計算

lambdas, eigenvectors = np.linalg.eig(Sigma)

固有値の大きい順にソートする

idx = lambdas.argsort()[::-1]
lambdas, eigenvectors = lambdas[idx], eigenvectors[idx]

変換する

transformed = X_ @ eigenvectors

これだけ

関数全体

def pca(X, n_components=2):
  """
  Args)
    X (np.array):
      データ行列
  Returns)
    transformed
  """
  # 共分散行列計算
  N = X.shape[0]
  X_ = X - X.mean(axis=0).reshape(1,-1)
  Sigma = 1/N * X_.T @ X_

  # 共分散行列の固有値、固有ベクトル計算
  ## lambdas[i]の固有ベクトルはeigenvectors[:,i]に対応する
  lambdas, eigenvectors = np.linalg.eig(Sigma)

  # 固有値が大きい順でソート＝射影先の分散が大きい順でソート
  idx = lambdas.argsort()[::-1]
  lambdas, eigenvectors = lambdas[idx], eigenvectors[idx]

  # 変換
  transformed = X_ @ eigenvectors

  return transformed[:,:n_components

実行例

2000次元の遺伝子発現量を特徴量としてもつ細胞のデータ(細胞数3102)を2次元に圧縮してプロットすると下図のようになった。何日目にサンプルされたかで色を分けて表示すると、最初の3日間と20日以降で異なるクラスタを形成していることが確認できる。