カシミアのニット

カシミヤのニット

日々のメモです

Selective Ensemble under Regularization Framework

Selective Ensemble under Regularization Framework

Li, N. and Zhou, Z.-H.: Selective Ensemble under Regularization Framework, Proc. MCS, pp. 293–303 (2009).

link.springer.com

実装はここに置いています.

github.com

アンサンブル枝刈り

アンサンブル枝刈りはensemble pruningとも呼ばれます.

アンサンブル枝刈りは訓練された複数の基本モデルが与えられた上で,その全てを結合するのではなく,部分集合を選択することです.

アンサンブル枝刈りにより,より小さなサイズのアンサンブルでより良い汎化性能を得ることが期待されます.

Tsoumakasらは,アンサンブル枝刈りを順序付けに基づく枝刈り,クラスタリングに基づく枝刈り,最適化に基づく枝刈りの3つのカテゴリーに分類しました*1

この詳しい話については後ほど別の記事にしようと思ってます.

今回はその中でも最適化に基づく枝刈りに属するRSE(regularized selective ensemble)についてです.

RSE(regularized selective ensemble)*2

RSEはLiとZhouによって提案されたアンサンブル枝刈りをQP問題へと帰着させる手法です.

この手法の良いところは,従来のアンサンブル枝刈り手法よりもアンサンブルに含まれる基本モデルの数が少ないにもかかわらず,その性能がいいという点です.

また,正則化項にグラフラプラシアン正則化を用いており,半教師あり学習の場面でも使用することができます.

それでは詳しく見ていきましょう.

{M}個のモデル{\{h_1, \ldots, h_M\}}に対し,アンサンブル結合重みベクトルを{\boldsymbol{w} = [w_1, \ldots, w_M]^{\mathrm{T}}}と定義します.

この時,{w_i \geq 0} かつ {\sum_{i=1}^{M}w_i = 1}です.

RSEは,以下の正則化リスク関数を最小化することにより{\boldsymbol{w}}を決定します.

\begin{align} R(\boldsymbol{w}) = \lambda V(\boldsymbol{w}) + \Omega(\boldsymbol{w}) \end{align}

ここで,{V(\boldsymbol{w})}は訓練データ{D = {(\boldsymbol{x}_1, y_1), \ldots, (\boldsymbol{x}_N, y_N)}}に対する誤分類の経験損失で, {\Omega(\boldsymbol{w})}正則化項であり, {\lambda}{V(\boldsymbol{w})}{\Omega(\boldsymbol{w})}の最小化における正則化パラメータを表す.

ヒンジ損失とグラフラプラシアン正則化項をそれぞれ経験損失と正則化として用いることにより,問題は下式で定式化されます.

\begin{align} & \underset{\boldsymbol{w}}{\text{minimize}} & & \boldsymbol{w}^{\mathrm{T}}{\boldsymbol{PLP}^{\mathrm{T}}}\boldsymbol{w} + \lambda\sum_{i=1}^N {\text{maximize}}(0, 1 - y_i\boldsymbol{p}_i^{\mathrm{T}}\boldsymbol{w})\\ & \text{subject to} & & \boldsymbol{1}^{\mathrm{T}}\boldsymbol{w} = 1, \; \boldsymbol{w} \geq \boldsymbol{0}. \label{eq:hinge_laplacian} \end{align}

ここで,{\boldsymbol{p}_i = (h_1(\boldsymbol{x}_i), \ldots, h_M(\boldsymbol{x}_i))^{\mathrm{T}}}は訓練データ{\boldsymbol{x}_i}に対する個々のモデルの予測を表し, {\boldsymbol{P} \in {{-1, +1}}^{M \times N}}は全訓練データに対する全モデルの予測を集めた予測行列で, {\boldsymbol{P}_{ij} = h_i(\boldsymbol{x}_j)}です.

{\boldsymbol{L}}は訓練データの近傍グラフ{G}の正規化グラフラプラシアンです.

グラフラプラシアンの詳細については,長くなってしまうので件の論文や他文献に譲ります.

行列{G}の重み付け隣接行列を{\boldsymbol{W}}で表し,{\boldsymbol{D}}{D_{ii} = \sum_{j=1}^N W_{ij}}の対角行列です.

この時,{\boldsymbol{L} = \boldsymbol{D}^{-\frac{1}{2}}(\boldsymbol{D} - \boldsymbol{W})\boldsymbol{D}^{-\frac{1}{2}}}となります.

上式の{maximize(\cdot)}は滑らかではので,スラック変数{{\boldsymbol{\xi}} = (\xi_1, \ldots, \xi_m)^{\mathrm{T}}}を導入することにより,上式は下式で書き表せます.

\begin{align} & \underset{\boldsymbol{w}}{\text{minimize}} & & \boldsymbol{w}^{\mathrm{T}}{\boldsymbol{PLP}^{\mathrm{T}}}\boldsymbol{w} + \lambda\boldsymbol{1}^{\mathrm{T}}\boldsymbol{\xi} \\ & \text{subject to} & & y_i\boldsymbol{p}_i^{\mathrm{T}}\boldsymbol{w} + \xi_i \leq 1, \; (\forall i = 1, \ldots, N) \\ & & & \boldsymbol{1}^{\mathrm{T}}\boldsymbol{w} = 1, \; \boldsymbol{w} \geq \boldsymbol{0}, \; \boldsymbol{\xi} \geq \boldsymbol{0}. \label{eq:hinge_laplacian_slack} \end{align}

この時,上式は標準的なQP問題となり,従来の最適化パッケージを用いて,効率的に解くことができます.

また,{\boldsymbol{1}^{\mathrm{T}}\boldsymbol{w} = 1, \boldsymbol{w} \geq \boldsymbol{0}}という制約は、スパース誘導性を持つ{l_1}ノルム制約となり,重み{\boldsymbol{w}}のいくつかの要素を強制的にゼロにします.

この特徴により,従来のアンサンブル枝刈り手法よりもアンサンブルに含まれる基本モデルの数が少ないにもかかわらず,それなりの性能を出すことができます.

また,簡単な例を代入すればわかりやすいのですが,正則化 {\boldsymbol{PLP}}の部分が似ているデータに対して違うラベルを割り当てると正則化するという振る舞いをします.

導出された結合重みベクトル{\boldsymbol{w}}を用いて,下式のようにアンサンブルの予測を得ます.

\begin{align} H(\boldsymbol{x}) = \sum_{w_i > 0}w_ih_i(\boldsymbol{x}) \quad\text{(RSE-w)} \label{eq:rse-w} \end{align}

また,下式のように{\boldsymbol{w}}の要素がゼロでない候補モデルの投票により予測を決定する選択的アンサンブルも提案されています.

\begin{align} H(\boldsymbol{x}) = \sum_{w_i > 0}h_i(\boldsymbol{x}) \qquad\text{(RSE)} \label{eq:rse} \end{align}

少なくともerrorに関してはなかなかいい結果です.

f:id:lapis_zero09:20180321163932p:plain

errorだけの比較なのは,あまりいい感じはしませんが,アンサンブルの目指す所的にいいのかな

実装

実装はここに置いています.

github.com

最適化を行う前に,グラフラプラシアンに用いる行列{\boldsymbol{W}}を導出する必要があります.

これは,pythonでは遅すぎたのでJavaで実装しています.

しかし,訓練データはnpz形式で保存していたので,pythonでデータを読み込んでpy4j*3で起動したJavaサーバで計算を行っています.

    /**
     * get the matrix W for Laplacian
     * @return w_link filename
     * @throws Exception
     */
    public String[] getKernelLinkMatrix(byte[] data, String fold) throws Exception {
        System.out.println("\t\t[*] Hello from Java!");
        // get the link matrix
        // create java matrix from numpy
        Instances d = createFromPy4j(data);
        Matrix LMat = new Matrix(m_numInst, m_numInst);
        RBFKernel krl = new RBFKernel();
        double g = 0.5 / m_numAttr;
        krl.setGamma(g);
        krl.buildKernel(d);

        System.out.println("\t\t[*] calculating w_link on Java...");
        for (int i = 0; i < m_numInst; i++) {
            for (int j = 0; j <= i; j++) {
                double vt = krl.eval(i, j, d.instance(i));
                LMat.set(i, j, vt);
                LMat.set(j, i, vt);
            }
        }
        d = null;
        krl = null;
        System.out.println("\t\t[-] calc done...");

        return returnFilename(LMat, m_numInst, fold);
}

ラプラシアン行列の計算や重みの最適化はMatlabで行いました.

本件ではpythonよりも速く,また,最適化パッケージであるmosek*4の扱いが簡単だったためです.

mosekは論文でも採用されており,QP問題に関しては最速らしいです.

function Lap = LaplacianMatrix(W)
% compute the laplacian matrix
%
% Input:
%       W: the line matrix
% Output:
%       Lap: the laplacian matrix


D = diag(sum(W,2));
N = D^(-0.5);
Lap = N *(D-W)* N;
function weight = compute_weight(lambda, Y, Prd, Q, dirname);
% compute the weights for the base classifiers
%
% Input:
%       Cfr: the ensemble classifier, which contains multiple base
%       classifiers
% Output:
%       weight: the weights for the base classifiers in the ensemble


% M : set the size of classifiers
% N : set the size of training data
[M, N] = size(Prd)
NumVar = M + N;

% ---- Prepare QP ----
%equal constaint

Aeq = [ones(1,M),zeros(1,NumVar - M)];
beq = 1;
% lower and upper bounds
lb = zeros(NumVar,1);
% Ax <= b
AP = Prd';
for id = 1:N
    AP(id, :) = AP(id, :) * Y(id);
end
AP = -1 .* AP;
A = [ AP, -1 * eye(N) ];
b = -1 * ones(N ,1);

% H
H = zeros(NumVar,NumVar);
H(1:M, 1:M) = Q;
% f
% lambda = 1.0;
f = lambda * [zeros(M,1);ones(N, 1)];

% Optimization Using QP - MOSEK
options = optimset('Display','off','TolFun', 1e-100);
x0 = quadprog(H,f,A,b,Aeq,beq,lb,[],[],options);
weight = x0(1:M);
weight((weight <= 1e-6)) = 0;
filename = sprintf('%sweight/weight_lambda_%d%s',dirname,lambda,'.csv')
disp(filename)
csvwrite(filename,weight);

簡単な実験

Phonemeデータセットに対して,5分割交差検証でRSEを使用しました.

SVM,KNN,Decision Tree,Bagged Decision Trees,Boosted Decision Trees,Boosted Decision Stumpsなどを様々なパラメータで767個のモデルをアンサンブルの候補としました.

このとき最もよかった候補モデルは,XGBoostではなくRandomForestで,mseが0.092937219730941698,f1が0.93484136670163342という値でした.

RSEは,mseが0.0859865470852018,f1が0.93981530924414192という結果だったので良くなっていますね.

best base model RSE
mse 0.0929372197309416 0.93484136670163342
f1 0.0859865470852018 0.93981530924414192

詳細は,

github.com

感想など

グラフラプラシアンを使っている以上うまく適合するデータとそうでないデータはもちろんあると思います.

また,この手法では,アンサンブル候補モデルに対する結合重みを二次計画問題と捉え最適解を導出しますが,候補モデルが一律に似たような性能を示している場合,二次計画問題となり得ず,解くことができません.

{l_1}ノルム制約を使って候補モデルの数をガンガン減らすのはとても理にかなっていると思ったので,回帰とかにも使えると嬉しいですね.

*1:Tsoumakas, G., Partalas, I. and Vlahavas, I.: An En- semble Pruning Primer, SUEMA, pp. 1–13 (2009).

*2:Li, N. and Zhou, Z.-H.: Selective Ensemble under Regularization Framework, Proc. MCS, pp. 293–303 (2009).

*3:https://www.py4j.org/

*4:https://www.mosek.com/