はじめに

機械学習コンペでもよく使われるXGBoostやLightGBMのバックグラウンドとなっている、決定木(Decision Tree)を用いた勾配ブースティング(Gradient Boosting)について、ヒストグラムを用いた手法を中心に、原論文を参考にしながら実装してみました。
勾配ブースティングの理論背景、決定木で用いるsplitterの決定基準などについても解説しました。
以下が参考にした論文です。

最急降下法

まずは勾配ブースティングの基本になっている、最急降下法の考え方の復習から始めたいと思います。
制約なし最適化問題の解法の一つである最急降下法は、ある関数(損失関数)の極小値を与える変数を、関数の傾きから求める、最も単純で手軽な最適化手法です。

L = l(\alpha)

という関数があったときに、これを現在の値よりも小さくするためには、$\alpha$を

\alpha \rightarrow \alpha-\gamma\frac{dl(\alpha)}{d\alpha}

と更新します。
これは、値$\alpha$において関数$l(\alpha)$が正の傾きを持つ時は、$\alpha$を負の方向に、逆に関数$l(\alpha)$が負の傾きを持つ時は、$\alpha$を正の方向に、変化させることで、$l(\alpha)$の値が小さくなるからです。変化させる値の大きさはパラメータ$\gamma$で調整します。

勾配ブースティング

その名の通り勾配ブースティングはブースティングの一種なので、複数の弱学習器を結合させることで強学習器を作ることを目的とします。
各予測器は逐次的に学習され、前の予測器の結果を随時修正していきます。
t番目の予測器を$f_t$と書くと、最終的な予測器$F_T$は

F_T = \sum_{t} \gamma_t f_t

と表されます。

（本記事では入力$x_i$から実数$y_i$を予測する回帰問題を扱います。）
最終的な目標は損失関数$L^T$を最小化することです。

L^T = \sum_{i} l(y_i, F_T(x_i))

$l$として今回は二乗誤差に係数$\frac{1}{2}$をかけたものを用いました。
t-1までの結果が与えられているとすると、その時の損失関数は

L^{t-1} = \sum_{i} l(y_i, F_{t-1}(x_i))

で与えられます。
ここで、この損失関数$L^{t-1}$を現在の値よりも小さくするためには、関数$F_{t-1}(x_i)$を変数と見て、各サンプルの損失関数$l$を変数$F_{t-1}(x_i)$に関して小さくすれば良いです。
そのためには、上記の最急降下法の考え方を用いて、

F_{t-1}(x_i) \rightarrow F_{t-1}(x_i) - \gamma \frac{\partial l}{\partial F_{t-1}(x_i)}

とすれば良いです。
勾配ブースティングでは逐次的に予測器を学習していくので、実際には矢印の左辺は$F_t(x_i)$になります。

F_{t}(x_i) = F_{t-1}(x_i) - \gamma \frac{\partial l}{\partial F_{t-1}(x_i)}

上式と$F_t(x_i)$の定義式

F_{t}(x_i) = F_{t-1}(x_i) + \gamma f_t(x_i)

を見比べてみると

f_t(x_i) = -\frac{\partial l}{\partial F_{t-1}(x_i)}

となり各予測器は、前のイテレーションまでの予測の総和の値 $F_{t-1}(x_i)$ における損失関数 $l$ の微分にマイナスをつけたものを目標値として学習させれば良いことが分かります。

r_i = -\frac{\partial l}{\partial F_{t-1}(x_i)}

とおくと、
各イテレーションにおいて学習させる予測器$f_t$の学習には次のコスト関数

C^t = \sum_{i}(r_i-f_t(x_i))^2

を最小化させれば良いということになります。個々のサンプルのコスト関数には二乗誤差を用います。
パラメータ$\gamma$は以下のように求めます。

\newcommand{\argmin}{\mathop{\rm argmin}\limits}

\gamma^* = \argmin_{\gamma} \sum_{i} l(y_i, F_{t-1}(x_i)+\gamma f_t(x_i))

t=1では$F_{t-1}$は無いので$r_i=y_i$と初期化します。

これで各イテレーションで予測器をどのように学習させれば良いのかが分かりました。
本記事では予測器$f_t$として決定木(Decision Tree)を用います。

決定木の学習

決定木はノードと二つの分岐を持つ二分木です。各ノードでは、いずれかの特徴量とその基準値が与えられ、各サンプルの対応する値がその基準値よりも小さければ左のノードへ、そうでなければ右のノードへ振り分けられます。

基準値(splitter)の候補の選び方

各ノードではある基準値(splitter)に基づいてサンプルを左右のノードに振り分けますが、その基準値はいくつかの候補の中から、ある指標によって選ばれます。その候補値の選び方はいろいろあるのですが、今回は代表的な二つを紹介します。

Exact Greedy
最も愚直で、精度は高いが時間がかかる手法です。一つの特徴量に着目すると、全ての訓練サンプル$x_i$をその特徴量の値でソートし、全てのサンプルの境界$s_i = (x_i+x_{i+1})/2 $をsplitterの候補とします。サンプル数を$N$とすると、$N-1$個の候補値ができることになります。この手法はscikit-learnやXGBoostで使われています。
Histogram Based
上記の手法よりも速く、ヒストグラムのビンやsplitterの数をうまく指定すれば精度もあまり劣らない手法です。$B$個のビンを持つヒストグラムがあったとき、ヒストグラム$H$は$H=\lbrace (p_1, m_1, r_1), ..., (p_B, m_B, r_B) \rbrace$ と表されます。$p_i$はビンの代表値（サンプルの平均値）、$m_i$はビンに含まれるサンプルの個数、$r_i$はビンに含まれるサンプルの$r$の総和です。このヒストグラムの値から、splitterの候補を計算します。着目する特徴量の範囲を$\tilde{B}$個のintervalに分割し、各intervalが大体$\frac{N}{\tilde{B}}$個のサンプルを含むようにします。このintervalの境界がsplitterの候補値になります。具体的な計算方法については下記で述べます。$\tilde{B}$は総サンプル数$N$に比べてかなり小さいため、splitterの候補数はExact Greedyに比べてかなり少なくなります。この手法はXGBoostやLightGBMで使われています。

Histogramの作成方法

Histogramの作成時には全てのサンプルを見る必要があり、Histogram Basedの手法において最も
時間を要する工程です。
まず、新たなサンプル$x_i$が与えられたとき、
ヒストグラムに$p_{B+1}=x_i$, $m_{B+1}=1$, $r_{B+1}=r_i$として新たに$(p_{B+1}, m_{B+1}, r_{B+1})$ を加えます。次に

i^*,j^* = \argmin_{i,j} |p_i-p_j|

を満たす$i, j$を求めて、その二つのビンをマージします。つまり、新たなサンプルも含めた全てのビンの中から、最も値の近い二つのビンを一つにまとめます。マージによってできる新たなビンは以下のようになります。

\displaylines{
p_{new} = \frac{m_ip_i+m_jp_j}{m_i+m_j} \\
m_{new} = m_i + m_j\\
r_{new} = r_i + r_j
}

Histogram作成の実装

Histogram作成のコードは以下になります。

    def _make_hist(self, x, r):        
        hist = []
        for i in range(len(x)):
            if i+1 <= self.num_bins:
                hist.append((x[i], 1, r[i]))
                if i+1==self.num_bins:
                    hist = sorted(hist, key=lambda tup:tup[0])                        
                    p_diffs = []
                    for j in range(len(hist)-1):
                        p_diffs.append(hist[j+1][0] - hist[j][0])
                    p_diffs = np.array(p_diffs, dtype=np.float64)
            else:
                if x[i] < hist[0][0]:
                    p_diffs = np.insert(p_diffs, 0, hist[0][0]-x[i])
                    hist.insert(0, (x[i], 1, r[i]))
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                    
                elif hist[0][0] <= x[i] and x[i] < hist[-1][0]:
                    j = 0
                    while hist[j][0] <= x[i]:
                        j += 1
                    hist.insert(j, (x[i], 1, r[i]))
                    p_diffs[j-1] = x[i] - hist[j-1][0]
                    p_diffs = np.insert(p_diffs, j, hist[j+1][0]-x[i])
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                    
                else:
                    hist.insert(len(hist), (x[i], 1, r[i]))
                    p_diffs = np.insert(p_diffs, self.num_bins-1, hist[-1][0]-hist[-2][0])
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                                       
        return hist

関数_make_hist()でヒストグラムを作成します。
引数xは与えられたサンプル$X$から特定の特徴量$X[:, k]$のみを取り出したものです。
与えられたビン数に到達するまでサンプルを追加していき、規定数に達した段階で、いったんヒストグラムをソートしています。その後上で述べた方法を用いてサンプルが加えられるたびにビンをマージしていきます。
マージのコードは以下になります。

    def _merge_hist(self, min_ind, hist, p_diffs):
        p1, p2 = hist[min_ind][0], hist[min_ind+1][0]
        m1, m2 = hist[min_ind][1], hist[min_ind+1][1]
        r1, r2 = hist[min_ind][2], hist[min_ind+1][2]
        
        hist[min_ind] = ((m1*p1 + m2*p2)/(m1+m2), m1+m2, r1+r2)
        hist.pop(min_ind+1)
        
        if min_ind == 0:
            p_diffs = np.delete(p_diffs, min_ind+1)
            p_diffs[min_ind] = float(hist[min_ind+1][0] - hist[min_ind][0])
        elif min_ind < self.num_bins-1:
            p_diffs = np.delete(p_diffs, min_ind+1)
            p_diffs[min_ind-1] = float(hist[min_ind][0] - hist[min_ind-1][0])
            p_diffs[min_ind] = float(hist[min_ind+1][0] - hist[min_ind][0])
        elif min_ind == self.num_bins-1:
            p_diffs[min_ind-1] = float(hist[min_ind][0] - hist[min_ind-1][0])
            p_diffs = np.delete(p_diffs, min_ind)
            
        assert len(hist) == self.num_bins and len(p_diffs) == self.num_bins-1
        
        return hist, p_diffs

splittterの候補の作成方法

関数_make_hist()で作成されたヒストグラムを
用いてsplitterの候補を作成します。
以下ヒストグラムはソートされているものとします。
候補値を$u_1, u_2, ..., u_{\tilde{B}-1}$とすると、$u_1$の左側、各$u_i$と$u_j$の間、それから$u_{\tilde{B}-1}$の右側にはそれぞれ大体$N/\tilde{B}$個のサンプルが含まれることになります。また各ビンの周りには$m$個のサンプルが均等に配置されているものと仮定します。さらにある値$x$よりも左側にあるサンプルの個数を$sum_l(x)$と表すことにします。$u_j$を求めたい時、まず$sum_l(u_j) = s_j = \frac{j}{\tilde{B}} N $を計算します。次に、$sum_l(p_i) < s_j < sum_l(p_{i+1})$を満たす$i$を見つけます。
ここで$sum_l(p_i)$は以下の式で求められます。

sum_l(p_i) = \sum_{j=1}^{i-1} m_j + \frac{m_i}{2}

$d = s_j - sum_l(p_i)$とし、
$a = m_{i+1} -m_i, b=2m_i, c=-2d$ とするとサンプルの個数が台形の面積に比例することを用いたtrapezoid methodを利用して$u_j$は

\displaylines{
u_j = p_i+(p_{i+1}-p_i)z \\
z = \frac{-b+\sqrt{b^2-4ac}}{2a} 
}

のようにして求めることができます。
trapezoid methodの詳細については論文[2 ] の2.1節をご覧ください。

splitterの候補の作成の実装

splitterの候補の作成に関わるコードは以下になります。

    def _make_sps(self, hist_, m_all):
        assert len(hist_) == self.num_bins+2
        b_ = self.num_sps+1
        cand_sps = []
        
        for j in range(self.num_sps):
            s = float((j+1)/b_) * m_all
            cand_s = self._uniform(s, hist_)
            cand_sps.append(cand_s)
        
        return cand_sps

    def _uniform(self, s, hist_):
        ### find j which satisfies sum_m(p_j) < s < sum_m(p_j+1)
        cand_sp = None
        j =1
        #now_m = self._sum_m(hist_[j][0], hist_)
        now_m = 0
        while now_m <= s:
            if now_m == s:
                cand_sp = hist_[j][0]
                return cand_sp
            if j == 1:
                now_m = hist_[j][1]/2
                
            j += 1
                
            #now_m = self._sum_m(hist_[j][0], hist_)
            ### for faster calculation 
            now_m = now_m + (hist_[j-1][1]+hist_[j][1])/2

        ###trapezoid method
        j -=1
        delta = 1e-10 ## to avoid a==0
        a = hist_[j+1][1] - hist_[j][1] + delta
        b = 2 * hist_[j][1]
        now_m = now_m - (hist_[j][1]+hist_[j+1][1]) / 2
        d = s - now_m
        #d = s - self._sum_m(hist_[j][0], hist_)
        c = -2*d
        z = (-b+np.sqrt(b**2 - 4*a*c))/ (2*a) # check if float or not 
        cand_sp = hist_[j][0] + (hist_[j+1][0]-hist_[j][0]) * z
        return cand_sp

hist_は_make_hist()で作成したhistの左端に$x$の最小値を右端に$x$の最大値を代表値に持つビンを加えたものです。これは求める$u_j$がhist中のビンの代表値の最小値よりも小さい、または最大値よりも大きい時のために必要な措置になります。
uniform()では $p_i<u_j<p_{i+1}$を満たす$i$(コード中では$j$)を求めたのち、上記のtrapezoid methodによって$u_j$を計算しています。

最適なsplitterの決定基準

Exact GreedyなりHistogram Basedなりでsplitterの候補を提案した後は、その中から最適なsplitterを選び出す必要があります。splitterの決定基準についてはさまざまなものがありますが、本記事では論文 [ 1 ]で提案されているものを採用します。
一つの予測器の学習に用いられるコスト関数$C$は下記のものでした。

C = \sum_{i} (r_i-f(x_i))^2

候補値の中からsplitterを一つ選び出し、それによってサンプルを左右に振り分けたとし、それぞれのインデックスの集合を$I_l^s$, $I_r^s$とすると,
上の$C$は

 C = \sum_{i\in I_l^s} (r_i-\bar{r_l^s})^2 + \sum_{i\in I_r^s} (r_i-\bar{r_r^s})^2

と書き換えられます。ここでは$\bar{r_l^s}$と$\bar{r_r^s}$は定数だということに注意してください。これは左右の子ノードを末端ノード(葉ノード)だとした時、決定木ではその末端ノードに達した全てのサンプルについての予測値は定数となるからです。
この定数値としては、二乗誤差を最小化する値である、全所属サンプルの平均値を取ることが多いです。
ここで
$\bar{r_l^s} = \frac{1}{|I_l^s|} \sum_{i\in I_l^s} r_i = \frac{l_s}{m_s} $とすると, $\bar{r_r^s} = \frac{l_{\infty}-l_s}{m_{\infty}-m_s}$
( $l_{\infty} = \sum_{i} r_i$, $m_{\infty} = N$ )
と表せられることに注意して、これらを上の$C$に代入して、式を展開し$s$に関わる部分のみを取り出すと、

C_s = -\frac{l_s^2}{m_s} - \frac{(l_{\infty}-l_s)^2}{m_{\infty}-m_s}

と変形できます。
$l_s$は$m_s=sum_l(s)$と同じ要領で今度は$r$について計算することで求められます。よって、このコスト関数を計算し、全ての候補をなめることで、コスト関数の値が最小となるsplitterを最適なsplitterとすることができます。
よって、最適なsplitterは以下のようにして求められることになります。

s^* = \argmin_s C_s

splitterを選定する関数の実装

一つの特徴量に着目して最適なsplitterを選び出す関数の実装は以下のようになります。

    def _fit_one_axis(self, x, r):
        # sps = splitters
        hist = self._make_hist(x, r)
                
        hist_ = hist.copy()
        min_ind = np.argmin(x)
        delta = 1e-10  ## for numerical calculation error 
        hist_.insert(0, (x[min_ind]-delta, 1, r[min_ind]))
        max_ind = np.argmax(x)
        hist_.append((x[max_ind]+delta, 1, r[max_ind]))
        
        cand_sps = self._make_sps(hist_, len(x))
        min_ind = None
        min_loss = float('inf')
        m_all = len(x)
        r_all = self._sum_r(hist_[-1][0], hist_)
        for i in range(len(cand_sps)):
            sp = cand_sps[i]
            if sp >= hist_[-1][0]:
                break
            m_s = self._sum_m(sp, hist_)
            r_s = self._sum_r(sp, hist_)
            loss = -r_s**2/m_s - (r_all-r_s)**2/(m_all-m_s)
            if loss < min_loss:
                min_loss = loss
                min_ind = i
        return min_loss, cand_sps[min_ind]

histとsplitterの候補cand_spsを作ったのち、
for文において、全ての候補について$m_s$, $r_s$(本文中では$l_s$)を計算し、コスト関数を求め、それが最小となるsplitterを返しています。ここで、$loss$の値まで返しているのは、今度は全ての特徴量の中からコストが最小となるものを選ぶ必要があるからです。

モデル全体の実装

以下がExact Greedy手法を用いたNodeGreedyモデルとHistogram Based手法を用いたNodeHistモデル(local)、NodeHistGlovalモデル(gloval)の実装になります。localとglovalの違いについては下記で説明しています。

import numpy as np 
import time
import copy

Exact Greedy

class NodeGreedy:   ##Exact Greedy method enumerate all feature and all splitter
    def __init__(self, d, max_depth, min_sample):
        self.d = d
        self.max_depth = max_depth
        self.min_sample = min_sample
        self.feature_ind = None
        self.splitter = None
        self.left_node = None
        self.right_node = None
        self.is_leaf = False
        self.r_pred = None
        
    def _fit_one_axis(self, x, r):
        
        #start = time.time()
        sorted_ind = np.argsort(x)
        #end = time.time()
        #print('   sorting takes {} seconds'.format(end-start))
        
        sorted_x = x[sorted_ind]
        sorted_r = r[sorted_ind]
        r_all = np.sum(sorted_r)
        m_all = len(x)
        losses, splits = [], []
        for i in range(len(sorted_x)-1):
            
            #start = time.time()
            
            split = (sorted_x[i]+sorted_x[i+1]) / 2
            ## number of samples after splitted should larger than min_sample
            if len(sorted_x[:i+1]) < self.min_sample or len(sorted_x[i+1:]) < self.min_sample:
                continue
            #if i+1 < self.min_sample or m_all- (i+1) < self.min_sample:
            #   continue
            
            ### deal with dummy feature
            #if np.sum(sorted_x==0) + np.sum(sorted_x==1) == len(sorted_x):
            #    if split!=0.5:
            #       continue
            
            r_s = np.sum(sorted_r[:i+1])
            m_s = i+1
            loss = -r_s**2 / m_s - (r_all-r_s)**2 / (m_all-m_s)
            losses.append(loss)
            splits.append(split)
        
        assert len(losses)!=0
        
        ind = np.argmin(losses)
        min_loss = losses[ind]
        best_split = splits[ind]
        return min_loss, best_split
            
        
    def fit(self, X, r):
        #print(' depth: {}, sample num: {}'.format(self.d, len(X)))
        if len(X) < 2*self.min_sample or self.d==self.max_depth:
            self.is_leaf = True
            self.r_pred = np.mean(r)
            return
        n_sample = X.shape[0]
        n_feature = X.shape[1]
        losses, splits = [], []
        for k in range(n_feature):
            #print(' {} th feature'.format(k+1))
            #start = time.time()
            loss, split = self._fit_one_axis(X[:,k], r)
            #end = time.time()
            #print('  takes {} time'.format(end-start))
            losses.append(loss)
            splits.append(split)
            
        assert len(losses) == n_feature and len(splits) == n_feature
        
        self.feature_ind = np.argmin(losses)
        self.splitter = splits[self.feature_ind]
        
        left_samples_ind = np.where(X[:, self.feature_ind] < self.splitter)[0]
        right_samples_ind = np.where(X[:, self.feature_ind] >= self.splitter)[0]
        
        """
        left_samples_ind = []
        right_samples_ind = []
        for i in range(n_sample):
            if X[i, self.feature_ind] < self.splitter:
                left_samples_ind.append(i)
            else:
                right_samples_ind.append(i)
        """
        
        self.left_node = NodeGreedy(self.d+1, self.max_depth, self.min_sample)
        self.right_node = NodeGreedy(self.d+1, self.max_depth, self.min_sample)
           
        # depth first         
        self.left_node.fit(X[left_samples_ind, :], r[left_samples_ind])
        self.right_node.fit(X[right_samples_ind, :], r[right_samples_ind])
        
        
    def predict(self, x):
        if self.is_leaf:
            return self.r_pred
        r_pred = None
        if x[self.feature_ind] < self.splitter:
            r_pred = self.left_node.predict(x)
        else:
            r_pred = self.right_node.predict(x)
        return r_pred

Histogram Based (local)

class NodeHist(NodeGreedy):  ## split supposing is local way
    def __init__(self, d, max_depth, min_sample, num_splitter=100, num_bins=100 ):
        super().__init__(d, max_depth, min_sample)
        self.num_sps = num_splitter
        self.num_bins = num_bins
        
    def _merge_hist(self, min_ind, hist, p_diffs):
        p1, p2 = hist[min_ind][0], hist[min_ind+1][0]
        m1, m2 = hist[min_ind][1], hist[min_ind+1][1]
        r1, r2 = hist[min_ind][2], hist[min_ind+1][2]
        
        hist[min_ind] = ((m1*p1 + m2*p2)/(m1+m2), m1+m2, r1+r2)
        hist.pop(min_ind+1)
        
        if min_ind == 0:
            p_diffs = np.delete(p_diffs, min_ind+1)
            p_diffs[min_ind] = float(hist[min_ind+1][0] - hist[min_ind][0])
        elif min_ind < self.num_bins-1:
            p_diffs = np.delete(p_diffs, min_ind+1)
            p_diffs[min_ind-1] = float(hist[min_ind][0] - hist[min_ind-1][0])
            p_diffs[min_ind] = float(hist[min_ind+1][0] - hist[min_ind][0])
        elif min_ind == self.num_bins-1:
            p_diffs[min_ind-1] = float(hist[min_ind][0] - hist[min_ind-1][0])
            p_diffs = np.delete(p_diffs, min_ind)
            
        assert len(hist) == self.num_bins and len(p_diffs) == self.num_bins-1
        
        return hist, p_diffs
        
    
    def _make_hist(self, x, r):
        hist = []
        for i in range(len(x)):
            if i+1 <= self.num_bins:
                hist.append((x[i], 1, r[i]))
                if i+1==self.num_bins:
                    hist = sorted(hist, key=lambda tup:tup[0])                        
                    p_diffs = []
                    for j in range(len(hist)-1):
                        p_diffs.append(hist[j+1][0] - hist[j][0])
                    p_diffs = np.array(p_diffs, dtype=np.float64)
            else:
                if x[i] < hist[0][0]:
                    p_diffs = np.insert(p_diffs, 0, hist[0][0]-x[i])
                    hist.insert(0, (x[i], 1, r[i]))
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                    
                elif hist[0][0] <= x[i] and x[i] < hist[-1][0]:
                    j = 0
                    while hist[j][0] <= x[i]:
                        j += 1
                    hist.insert(j, (x[i], 1, r[i]))
                    p_diffs[j-1] = x[i] - hist[j-1][0]
                    p_diffs = np.insert(p_diffs, j, hist[j+1][0]-x[i])
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                    
                else:
                    hist.insert(len(hist), (x[i], 1, r[i]))
                    p_diffs = np.insert(p_diffs, self.num_bins-1, hist[-1][0]-hist[-2][0])
                    min_ind = np.argmin(p_diffs)
                    hist, p_diffs = self._merge_hist(min_ind, hist, p_diffs)
                                       
        return hist
    
    def _calc_a(self, u, p_i, p_i_1, x_i, x_i_1):
        x_u = x_i + ((x_i_1-x_i)/(p_i_1-p_i))*(u-p_i)
        a = ((x_i+x_u)*(u-p_i))/2
        return a
    
    def _sum_r(self, u, hist_):
        if u == hist_[-1][0]:
            ret = 0
            for j in range(1, len(hist_)-1):
                ret += hist_[j][2]
            return ret
        
        
        j = 0
        while hist_[j][0] <= u:
            j += 1
        j -= 1
        # p_j <= u < p_j+1
        a_u = self._calc_a(u, hist_[j][0], hist_[j+1][0], hist_[j][2], hist_[j+1][2])
        a_p_i_1 = self._calc_a(hist_[j+1][0], hist_[j][0], hist_[j+1][0], hist_[j][2], hist_[j+1][2])
        
        
        if j==0:
            R = hist_[j+1][2]/2
            return a_u / a_p_i_1 * R
        elif j==len(hist_)-2:
            R = hist_[j][2]/2
            r_l = 0
            for k in range(1, len(hist_)-2):
                r_l += hist_[k][2]
            r_l += hist_[j][2]/2
            return r_l + a_u/a_p_i_1 * R
        else:
            R = ( hist_[j][2]+hist_[j+1][2]) / 2
            r_l = 0
            for k in range(1, j):
                r_l += hist_[k][2]
            r_l += hist_[j][2]/2
            return r_l + a_u/a_p_i_1 * R
        
        
    
    def _sum_m(self, u, hist_):
        if u == hist_[0][0]:
            return 0
        elif u == hist_[-1][0]:
            ret = 0
            for j in range(1, len(hist_)-1):
                ret += hist_[j][1]
            return ret
        
        j = 0
        while hist_[j][0] <= u:
            j += 1     
        
        j -= 1
        # p_j <= u < p_j+1
        a_u = self._calc_a(u, hist_[j][0], hist_[j+1][0], hist_[j][1], hist_[j+1][1])
        a_p_i_1 = self._calc_a(hist_[j+1][0], hist_[j][0], hist_[j+1][0], hist_[j][1], hist_[j+1][1])
        
        
        if j==0:
            #M = hist_[j][1] + hist_[j+1][1]/2
            M = hist_[j+1][1]/2
            return a_u / a_p_i_1 * M
        elif j==len(hist_)-2:
            #M = hist_[j][1]/2 + hist_[j+1][1]
            M =  hist_[j][1]/2
            m_l = 0
            for k in range(1, len(hist_)-2):
                m_l += hist_[k][1]
            m_l += hist_[j][1]/2
            return m_l + a_u/a_p_i_1 * M
        else:
            M = ( hist_[j][1]+hist_[j+1][1]) / 2
            m_l = 0
            for k in range(1, j):
                m_l += hist_[k][1]
            m_l += hist_[j][1]/2
            return m_l + a_u/a_p_i_1 * M
    
    
    def _uniform(self, s, hist_):
        ### find j which is sum_m(p_j) < s < sum_m(p_j+1)
        cand_sp = None
        j =1
        #now_m = self._sum_m(hist_[j][0], hist_)
        now_m = 0
        while now_m <= s:
            if now_m == s:
                cand_sp = hist_[j][0]
                return cand_sp
            if j == 1:
                now_m = hist_[j][1]/2
                
            j += 1
                
            #now_m = self._sum_m(hist_[j][0], hist_)
            # for faster calculation 
            now_m = now_m + (hist_[j-1][1]+hist_[j][1])/2
            
        j -=1
        delta = 1e-10 ## to avoid a==0
        a = hist_[j+1][1] - hist_[j][1] + delta
        b = 2 * hist_[j][1]
        now_m = now_m - (hist_[j][1]+hist_[j+1][1]) / 2
        d = s - now_m
        #d = s - self._sum_m(hist_[j][0], hist_)
        c = -2*d
        z = (-b+np.sqrt(b**2 - 4*a*c))/ (2*a) # check if float or not 
        cand_sp = hist_[j][0] + (hist_[j+1][0]-hist_[j][0]) * z
        return cand_sp
            
    
    def _make_sps(self, hist_, m_all):
        assert len(hist_) == self.num_bins+2
        b_ = self.num_sps+1
        cand_sps = []
        
        for j in range(self.num_sps):
            s = float((j+1)/b_) * m_all
            cand_s = self._uniform(s, hist_)
            cand_sps.append(cand_s)
        
        return cand_sps
        
                         
        
    # override _fit_one_axis of NodeGreedy
    def _fit_one_axis(self, x, r):
        # sps = splitters
        hist = self._make_hist(x, r)
                
        hist_ = hist.copy()
        min_ind = np.argmin(x)
        delta = 1e-10  ## for numerical calculation error 
        hist_.insert(0, (x[min_ind]-delta, 1, r[min_ind]))
        max_ind = np.argmax(x)
        hist_.append((x[max_ind]+delta, 1, r[max_ind]))
        
        cand_sps = self._make_sps(hist_, len(x))
        min_ind = None
        min_loss = float('inf')
        m_all = len(x)
        r_all = self._sum_r(hist_[-1][0], hist_)
        for i in range(len(cand_sps)):
            sp = cand_sps[i]
            if sp >= hist_[-1][0]:
                break
            m_s = self._sum_m(sp, hist_)
            r_s = self._sum_r(sp, hist_)
            loss = -r_s**2/m_s - (r_all-r_s)**2/(m_all-m_s)
            if loss < min_loss:
                min_loss = loss
                min_ind = i
        return min_loss, cand_sps[min_ind]
    
    def fit(self, X, r):
        print(' depth: {}, sample num: {}'.format(self.d, len(X)))
        if len(X) < 2*self.min_sample or self.d==self.max_depth or len(X)<self.num_bins:
            self.is_leaf = True
            self.r_pred = np.mean(r)
            return
        n_sample = X.shape[0]
        n_feature = X.shape[1]
        losses, splits = [], []
        for k in range(n_feature):
            #start = time.time()
            loss, split = self._fit_one_axis(X[:,k], r)
            #end = time.time()
            #print('  takes {} time'.format(end-start))
            losses.append(loss)
            splits.append(split)
            
        assert len(losses) == n_feature and len(splits) == n_feature
        
        self.feature_ind = np.argmin(losses)
        self.splitter = splits[self.feature_ind]
        
        left_samples_ind = np.where(X[:, self.feature_ind] < self.splitter)[0]
        right_samples_ind = np.where(X[:, self.feature_ind] >= self.splitter)[0]
                
        self.left_node = NodeHist(self.d+1, self.max_depth, self.min_sample, self.num_sps, self.num_bins)
        self.right_node = NodeHist(self.d+1, self.max_depth, self.min_sample, self.num_sps, self.num_bins)
           
        # depth first         
        self.left_node.fit(X[left_samples_ind, :], r[left_samples_ind])
        self.right_node.fit(X[right_samples_ind, :], r[right_samples_ind])

Histogram Based (gloval)
これまで(localな手法)は各ノードごとにsplitterの候補値を計算していたのですが、glovalな手法では深さ0で始めに候補値を提案した後、その後のノードでも同じ候補値（最適なsplitterで分割済み）を用います。

class NodeHistGloval(NodeHist):
    def __init__(self, d, max_depth, min_sample, num_splitter=100, num_bins=100, sps_dic={}):
        super().__init__(d, max_depth, min_sample, num_splitter, num_bins)
        self.sps_dic = sps_dic
             
    def _fit_one_axis(self, x, r, f_ind):
        # sps = splitters
        hist = self._make_hist(x, r)
        hist_ = hist
        min_ind = np.argmin(x)
        delta = 1e-10
        hist_.insert(0, (x[min_ind]-delta, 1, r[min_ind]))
        max_ind = np.argmax(x)
        hist_.append((x[max_ind]+delta, 1, r[max_ind]))
        if self.d ==0 :
            cand_sps = self._make_sps(hist_, len(x))                
            self.sps_dic[f_ind] = cand_sps.copy()
        else:
            cand_sps = self.sps_dic[f_ind].copy()
                        
            if len(cand_sps)==0:
                return float('inf'), None, None
            
            ind = len(cand_sps)-1
            while cand_sps[ind] > np.max(x):
                cand_sps = cand_sps[:ind]
                ind -= 1
                
            ind = 0
            while cand_sps[ind] < np.min(x):
                cand_sps = cand_sps[ind+1:]
        
            self.sps_dic[f_ind] = cand_sps.copy()
              
        min_ind = None
        min_loss = float('inf')
        m_all = len(x)
        r_all = self._sum_r(hist_[-1][0], hist_)
    
        for i in range(len(cand_sps)):
            sp = cand_sps[i]
            m_s = self._sum_m(sp, hist_)
            r_s = self._sum_r(sp, hist_)
            loss = -r_s**2/m_s - (r_all-r_s)**2/(m_all-m_s)
            if loss < min_loss:
                min_loss = loss
                min_ind = i
        return min_loss, cand_sps[min_ind], min_ind
    
    
    def fit(self, X, r):
        print(' depth: {}, sample num: {}'.format(self.d, len(X)))
        if len(X) < 2*self.min_sample or self.d==self.max_depth or len(X)<self.num_bins:
            self.is_leaf = True
            self.r_pred = np.mean(r)
            return
        n_sample = X.shape[0]
        n_feature = X.shape[1]
        losses, splits, min_inds_sp = [], [], []
        for k in range(n_feature):
            #print(' {} th feature'.format(k+1))
            #start = time.time()
            loss, split, min_ind_sp = self._fit_one_axis(X[:,k], r, k)
            #end = time.time()
            #print('  takes {} time'.format(end-start))
            losses.append(loss)
            splits.append(split)
            min_inds_sp.append(min_ind_sp)
            
        assert len(losses) == n_feature and len(splits) == n_feature
        
        if np.min(losses) == float('inf'):
            self.is_leaf = True
            self.r_pred = np.mean(r)
            return 
        
        self.feature_ind = np.argmin(losses)
        self.splitter = splits[self.feature_ind]
        min_ind_sp = min_inds_sp[self.feature_ind]
        
        
        left_samples_ind = np.where(X[:, self.feature_ind] < self.splitter)[0]
        right_samples_ind = np.where(X[:, self.feature_ind] >= self.splitter)[0]
        
                
        left_sps_dic = self.sps_dic.copy()
        left_sps_dic[self.feature_ind] = self.sps_dic[self.feature_ind][:min_ind_sp]
        right_sps_dic = self.sps_dic.copy()
        right_sps_dic [self.feature_ind] = self.sps_dic[self.feature_ind][min_ind_sp+1:]
                
        self.left_node = NodeHistGloval(self.d+1, self.max_depth, self.min_sample, self.num_sps, self.num_bins, left_sps_dic)
        self.right_node = NodeHistGloval(self.d+1, self.max_depth, self.min_sample, self.num_sps, self.num_bins, right_sps_dic)
                
        # depth first         
        self.left_node.fit(X[left_samples_ind, :], r[left_samples_ind])
        self.right_node.fit(X[right_samples_ind, :], r[right_samples_ind])

結果比較

sckit-learnで提供されているトイデータセット(California-housing)を用いて各手法の訓練時間、予測時間、回帰の精度を表すMSEの比較をしてみました。データセットの訓練サンプル数は$n=20640$, featureの数は8です。予測は同じ訓練サンプルを用いて行なっています。決定木の最大の深さは５、Histogram Basedの手法についてはsplitterの候補数は100, ビンの数は100としました。

n_iteration=1のとき

method	fit time	pred time	MSE
sckit-learn	0.091	0.002	0.4906
Exact Greedy	11.22	0.093	0.4911
Hist Based (local)	64.66	0.16	0.4994
Hist Based (gloval)	66.43	0.16	0.4997

scikit-learnの勾配ブースティングアルゴリズムはCythonで書かれているので、pythonで書いた自作のコードよりもかなり速かったです。精度的にはやはりExact Greedy手法を用いたscikit-learnとExact Greedyが良いことが分かります。この二つのモデルの精度の違いは、splitterの決定基準となるコスト関数の違いだと思われます。Exact Greedyは$O(n\log{n})$, Hist Basedは$O(n\log{B})$となるので、Hist Basedのモデルの方が通常速いのですが、今回Exact Greedyではソートにnumpyを用いているのに対して、Hist Basedではヒストグラムの作成時にfor文を用いて全てのサンプルを舐めているので、Hist Basedのモデルの方が遅いという結果になってしまいました。

n_iteration=10のとき

method	fit time	pred time	MSE
sckit-learn	0.836	0.0057	0.2339
Exact Greedy	113.12	0.942	0.2429
Hist Based (local)	677.02	1.00	0.2509
Hist Based (gloval)	665.98	1.04	0.2525

予測器の数を10まで増やすとブースティングの効果が表れ、どのモデルでも精度はかなり改善しています。自作のモデルたちはfittingに時間がかかり過ぎてしまっています。

最後に

最後に参考までにnumpyを用いたソートと単純なfor文とでそれぞれどれくらい時間がかかるのかを調べてみました。結果は以下のようになりました。($n=10^8$)

ソート・・・2.49s
for文・・・10.68s
コードは以下になります。

import numpy as np
import time
n=100000000
l = np.random.randint(1,10, size=n)
start = time.time()
l_sort = np.sort(l)
print('sort time:{} seconds'.format(time.time()-start))

start = time.time()
s = 0
for i in range(n):
    s += 1
print('for time:{} seconds'.format(time.time()-start))

お疲れ様でした。
長文にも関わらず、ここまで目を通していただきありがとうございました。全体を通して何らかの不備や、コードでもっとこうしたら速度が改善するなど、アドバイス等ありましたら教えていただけると助かります。