横須賀の某Prisonで働く研究者?のブログ

NAISTを出て横須賀の某社の研究所で研究者をしています.本ブログの内容は業務内容や所属企業とは一切関係ありません.

【論文解説/IEEE】Map-Reduceに対応した超高次元非線形特徴量選択/Ultra High-Dimensional Nonlinear Feature Selection for Big Biological Data

こんにちは,接点QBです. 今回は論文解説ネタです.

今回紹介する論文は Yamada et al. (2018) による Ultra High-Dimensional Nonlinear Feature Selection for Big Biological Data です. 本論文で提案されている手法であるLANDは,大雑把に言いうと,HSIC Lasso (Yamada et al. (2014))の高速版です. つまり,高次元小標本という設定で高速に特徴量選択をする手法です. 正規化されたHSIC (NHSIC) を計算する際にNystrom近似 (Scholkopf & Smola (2002)) を用いて, NHSICを構成する行列を低ランク近似で表現します. 低ランク近似された行列は並列に計算が可能なので,Map-Reduceで分散処理すると高次元のデータでも高速に処理できますよっていう話です.

問題設定

入力データを以下のように表します: \begin{equation} \mathbf{X}=\begin{bmatrix} \boldsymbol{x}_1 & \cdots & \boldsymbol{x}_{n} \end{bmatrix} =\begin{bmatrix} \boldsymbol{u}_{1} \\ \vdots \\ \boldsymbol{u}_d \end{bmatrix}^T \in \mathbb{R}^{d\times n}, \quad \boldsymbol{x}_i = \begin{bmatrix} x_{i1} \\ \vdots \\ x_{id} \end{bmatrix}, \quad \boldsymbol{u}_j = \begin{bmatrix} x_{1j} \\ \vdots \\ x_{dj} \end{bmatrix}^T. \end{equation} ここで,$\boldsymbol{x}_i$は第$i$サンプルを表し,$\boldsymbol{u}_{j}$はすべてのサンプルに渡る第$j$特徴量のベクトルを表しています. また,第$i$サンプルに対応する目的変数を$y_{i} \in \mathcal{Y}$ で表し, $\boldsymbol{y} = \begin{bmatrix} y_1 & \cdots & y_{n}\end{bmatrix}^T$ とします.

教師あり特徴量選択の目的は,$y$の予測に関係する$m << d$個の特徴量を発見することです. そのためにはいくつかの方法がありますが,Yamada et al. (2018) では HSIC Lasso with Least Angle Regression (LAND) という手法が提案されています. LANDはカーネル法を用いた手法なので,カーネル行列に対する表記を定義しておきます: \begin{align} \left[ \mathbf{K}^{(k)}\right]_{ij} & = K\left( u_{ki}, u_{kj}\right), \quad i,j = 1, \cdots, n \\ \left[ \mathbf{L}\right]_{ij} & = L \left( y_{i}, y_{j}\right), \quad i,j = 1, \cdots, n \end{align} ここで$K$と$L$はカーネル関数です. Yamada et al. (2018) ではカーネル関数としてガウシアンカーネルを用いています. なお,分類問題のときは \begin{equation} L(y, y') = \begin{cases} 1 / n_{y} & y = y' \\ 0 & \text{otherwise} \end{cases} \end{equation} を用いると良いです.なお,$n_{y}$は目的変数のクラス数です.

LAND

本題のLANDの話に入る前にHSIC Lasso (Yamada et al. (2014)) に触れておきます.

HSIC Lasso

HSIC Lasso は以下の最適化問題によって定式化されます: \begin{equation} \min_{\boldsymbol{\alpha}\in \mathbb{R}^d} \left\| \tilde{\mathbf{L}} - \sum_{k=1}^{d} \alpha_k \tilde{\mathbf{K}}^{(k)} \right\| ^2_F + \lambda \left\| \boldsymbol{\alpha}\right\| _1, \quad \text{s.t.} \quad \alpha_1, \cdots, \alpha_d \geq 0. \label{optimization problem hsiclasso} \end{equation} ここで,$\lambda \geq 0$ は正則化パラメータで,$\left\|\, \cdot\, \right\|_1$は$\ell_1$ノルム, $\left\|\, \cdot \, \right\|_F$はフロベニウスノルムです. また,$\tilde{\mathbf{K}}, \tilde{\mathbf{L}}$は以下を満たすカーネル行列です: \begin{align} & \boldsymbol{1}_n^T \tilde{\mathbf{K}} \boldsymbol{1}_n = \boldsymbol{1}_n^T \tilde{\mathbf{L}} \boldsymbol{1}_n^T = 0, \\ & \left\| \tilde{\mathbf{K}}\right\|_{F}^2 = \left\| \tilde{\mathbf{L}}\right\|_{F}^2 = 1. \end{align}

LAND

LANDはEq. (\ref{optimization problem hsiclasso}) をleast angle regressionで解く手法の名前です. LANDとHSIC Lassoの違いはEq. (\ref{optimization problem hsiclasso}) の解き方だけです.

Eq. (\ref{optimization problem hsiclasso}) の解により,最もアウトプットとの関連性が高く, 余剰が最も少なくなるような特徴量を選択することが可能になります. その説明のために,Eq. (\ref{optimization problem hsiclasso}) の目的関数を以下のように変形します: \begin{equation} 1 - 2\sum_{k=1}^{d} \alpha_{k} \textrm{NHSIC}(\boldsymbol{u}_k, \boldsymbol{y}) + \sum_{k=1}^{d}\sum_{l=1}^{d} \alpha_{k} \alpha_{l} \textrm{NHSIC}(\boldsymbol{u}_k, \boldsymbol{u}_l). \label{optimization problem2 hsiclasso} \end{equation} ここで,$\textrm{HSIC}(\boldsymbol{u}, \boldsymbol{y}) = {\rm tr} \left( \tilde{\mathbf{K}}\tilde{\mathbf{L}}\right)$はHSICのnormalized versionです. Eq. (\ref{optimization problem2 hsiclasso}) が何を言っているのかは,以下のように場合を分けて考えると分かりやすいと思います.

  1. $\boldsymbol{y}$が第$k$特徴量$\boldsymbol{u}_k$に依存する場合:$\textrm{NHSIC}(\boldsymbol{u}_k, \boldsymbol{y})$は大きくなり, したがって(Eq. (\ref{optimization problem2 hsiclasso)を最小化するので)$\alpha_k$も大きくなります. つまり,第$k$特徴量は関連性の高い特徴量として選択されることになります.
    1. $\boldsymbol{u}_k$に$\boldsymbol{u}_l$が依存する場合:Eq. (\ref{optimization problem2 hsiclasso}) 第3項のNHSICが大きくなるので, 全体を最小化するために$\alpha_l$が小さくなります. つまり,余剰な第$l$特徴量は選択されなくなります.
    2. $\boldsymbol{u}_k$に$\boldsymbol{u}_l$が依存しない場合: Eq. (\ref{optimization problem2 hsiclasso}) 第3項のNHSICは$0$に近い値となり,最適化問題への影響が小さくなります. また,この場合はEq. (\ref{optimization problem2 hsiclasso}) 第2項で$k=l$のケースによって第$l$特徴量が選択されるべきかを考えることになります.
  2. $\boldsymbol{y}$が第$k$特徴量$\boldsymbol{u}_k$と独立の場合:$\textrm{NHSIC}(\boldsymbol{u}_k, \boldsymbol{y})$は$0$に近くなり, したがって$\alpha_k$も小さくなります. つまり,第$k$特徴量は目的変数に対して関連性の低い特徴量として選択されないことになります. また,この場合は$\alpha_k$が$0$に近い値になるので,第3項の最適化問題全体への影響は小さくなります.

LANDは実際は反復的に目的変数に関連性の高い特徴量選択を行いますが, そのためには関連性の高い特徴量集合を作成し,その中から余剰な特徴量を取り除くという2つの段階が必要になります. このことから,特徴量の選択スコアを以下のように定義します: \begin{equation} c_k = \textrm{NHSIC}\left( \boldsymbol{u}_{k}, \boldsymbol{y}\right) - \sum_{i: \alpha_{i}>0}\alpha_{i} \textrm{NHSIC}\left( \boldsymbol{u}_k, \boldsymbol{u}_i \right) \end{equation}

ここまでは少し理論的な話をしてきましたが,Eq. (\ref{optimization problem hsiclasso}) を解くには計算コストが結構かかります. メモリだとカーネル行列$\tilde{\mathbf{K}}^{(k)}$をおいておくために$O(dn^2)$のメモリが必要になりますし, 計算時間的な観点から見ると,$O(mdn^3)$必要になります. この辺は「カーネル法あるある」なんですが,Yamada et al. (2018) ではこの問題の解決に取り組んでいます. 具体的にはNystrom approximationを使って,カーネル行列を行列の低ランク近似を用いて表します: \begin{align} \textrm{NHSIC}(\boldsymbol{y}, \boldsymbol{y}) &= {\rm tr}\left( \tilde{\mathbf{K}}\tilde{\mathbf{L}}\right) \approx {\rm tr} \left( \mathbf{F}\mathbf{F}^T \mathbf{G}\mathbf{G}^T \right), \\ \mathbf{F} & = \frac{ \mathbf{\Gamma}\mathbf{K}_{nb}\mathbf{K}_{bb}^{-1/2} }{ \left\{ {\rm tr}\left[ \left( \mathbf{K}^{-1/2}_{bb}\mathbf{K}_{nb}^{T}\mathbf{K}_{nb}\mathbf{K}_{bb}^{-1/2}\right)^{2} \right] \right\}^{1/4} }, \\ \mathbf{G} &= \frac{ \mathbf{\Gamma}\mathbf{L}_{nb}\mathbf{L}_{bb}^{-1/2} }{ \left\{ {\rm tr}\left[ \left( \mathbf{L}^{-1/2}_{bb}\mathbf{L}_{nb}^{T}\mathbf{L}_{nb}\mathbf{L}_{bb}^{-1/2}\right)^{2} \right] \right\}^{1/4} } \end{align} ここで,$\mathbf{F}\mathbf{F}^T$は$\tilde{\mathbf{K}}$の低ランク近似で, $\mathbf{K}_{nb}\in \mathbb{R}^{n\times b},\, \left[ \mathbf{K}_{nb}\right]_{ij} = K\left( u_{i}, u_{b,j}\right)$ を満たします. ただし,ここでの$\boldsymbol{u}_{b} \in \mathbb{R}^b$はバイアスベクトルです. 他の行列についても同様です. また,$b$は$\mathbf{K}_{nb}$と$\mathbf{L}_{nb}$の上界です. 一般的に,$b$を大きくするほど近似の精度は上がりますが,計算コストやメモリコストがかかります. さらに,各特徴量に対する$\left\{ \mathbf{F}_{k}\right\}_{k=1}^{d}$は並列計算が可能です. そのため,高次元データに対してはMap-Reduceを用いた並列処理が威力を発揮します.

まとめ

本記事で紹介した論文は,IEEE TRANSACTIONS OF KNOWLEDGE AND DATA ENGINEERING, VOL. 30, NO. 7に記載されているので, あまりすべてを勝手に無料で紹介してしまうのはマズいと思ったので,今回は手法の大枠の説明だけにとどめました. 原著にはアルゴリズム擬似コードや(当然ながら)既存手法との比較結果等が記載されていますので, 本記事を読んで興味を持たれた方は是非そちらを御覧ください.