DecisionTreeClassifierの中を見てみる

概要

機械学習をするためのメソッドは数多くありますが、kaggleなどで人気が高いのはXGBMやLightGBMなどの決定木を利用するものです。今回はそんな決定木が学習をするアルゴリズムを調べてみました。

詳細の前に

決定木とはなんぞや、についてはこの記事などを読んでみてください。

また、以下の前提を設けます。

・scikit-learnのDecisionTreeClassifierのfit()メソッドのソースコード を調べてみる。

・予測する変数は1種類

・サンプル毎の重みはつけない

・不純度の指標はエントロピー

全体像

1 stackに学習データを入れる

2 stackの先頭のデータ（X）を取り出す。

3 Xの特徴量（fi）を1つ選ぶ

4 Xを区切る位置（pos）を選ぶ

5 そのように区切った場合の不純度の減少量を求める

6 区切り位置をずらす

7 4に戻る

8 特徴量を変える

9 3に戻る

10 不純度の減少量が最大となる特徴量、区切り位置でXを分割

11 これ以上分割できない or 不純度が減少しなくなったら終了。

　 そうでない場合は分割した各々をstackに入れる

12 2に戻る

13 stackが空になったら終了

不純度というのはデータがどれだけ綺麗に分かれているかの指標で、小さい方が綺麗に分かれているということになります。不純度の計算方法はいくつかありますが、今回はエントロピーを例にしたいと思います。（エントロピーについては上記記事を参照）

全体像をコードで

BaseDecisionTreeのfit()メソッドを見てみます。

最初の方はパラメータのチェックで、実質TreeBuilderのbuild()メソッドを呼んでいるだけです。

sklearn/tree/_classes.py

class BaseDecisionTree(MultiOutputMixin, BaseEstimator, metaclass=ABCMeta):

def fit(self, X, y, sample_weight=None, check_input=True,X_idx_sorted="deprecated"):

//略

builder.build(self.tree_, X, y, sample_weight)

TreeBuilderのbuild()メソッドで上記の処理が行われます。

ここでは、DepthFirstTreeBuilderクラスのbuild()メソッドを見てみます。

sklearn/tree/_tree.pyx

cdef class DepthFirstTreeBuilder(TreeBuilder):
"""Build a decision tree in depth-first fashion."""

cpdef build(self, Tree tree, object X, np.ndarray y,                np.ndarray sample_weight=None):
"""Build a decision tree from the training set (X, y)."""

//略　変数の定義

with nogil:
# push root node onto stack  
    rc = stack.push(0, n_node_samples, 0, _TREE_UNDEFINED, 0, INFINITY, 0)
    //1 stackに学習データを入れる
    if rc ==-1:
    # got return code -1 - out-of-memory
    with gil:
        raiseMemoryError()
    while not stack.is_empty():                    
        stack.pop(&stack_record) //2 stackの先頭データを取り出す
        
        //略　データを取り出す処理

先頭でstackにデータ全体をpushし、stackにデータがある間はそのデータを取り出し、処理するというループが走ります。

        //略　データを取り出す処理
        is_leaf = (is_leaf or 
                  (impurity <= min_impurity_split))
                  
        // 3-10 データを分割
        if not is_leaf:                                
            splitter.node_split(impurity, &split, 　　　　　　　　　　　　　　　　　　　&n_constant_features) 
        # If EPSILON=0 in the below comparison, float precision
        # issues stop splitting, producing trees that are
        # dissimilar to v0.18                    
            is_leaf = (is_leaf or split.pos >= end or                                
                      (split.improvement + EPSILON < 
                       min_impurity_decrease)) 
            //これ以上分割できるかチェック
        node_id = tree._add_node(parent, is_left, is_leaf,   
                                 split.feature, split.threshold, 
                                 impurity, n_node_samples,                                         
                                 weighted_n_node_samples)

splitter.node_splitが上記3-10のデータを分割する特徴量と位置を探して最適なところで分割する処理です（詳細は後述します）。

tree.add_nodeで決定木を更新します。

#略
if not is_leaf:
    # Push right child on stack
    rc = stack.push(split.pos, end, depth +1, node_id, 0,                                    
                    split.impurity_right, n_constant_features)
    if rc ==-1:
        break
    # Push left child on stack   
    rc = stack.push(start, split.pos, depth +1, node_id, 1,                                    
                    split.impurity_left, n_constant_features)
    if rc ==-1:
        break
    

1個上のコードでis_leafの確認を行います。より良い分割位置が見つからなかった（split.pos >= end）もしくは不純度の減少量が小さい

（split.improvement + EPSILON < min_impurity_decrease）場合はstackには何も追加しません。そうでない場合は分割した2つのデータをそれぞれstackに追加します（11）。

これでwhileループは終了です。is_leafでない場合は分割したデータそれぞれについて同じ手順で分割を繰り返していきます。

データの分割

それではnode_splitメソッドの処理を見てみようと思います。分割を行うsplitterクラスも複数あるのですが、今回は最適な分割を探すBestSplitterのnode_splitメソッドを見てみたいと思います。

sklearn/tree/_splitter.pyx

cdef class BestSplitter(BaseDenseSplitter):
"""Splitter for finding the best split."""

本メソッドでは「全体像」のところで説明したように、データを分割するのに最適な特徴量と分割位置を探します。

特徴量のループ

まずは特徴量についてのループを回しますが、ここでは計算量節約のために1つずつ順番に調べるのではありません。与えられたデータサンプルで、データの値がほとんど変わらない場合はその特徴量については分割位置の探索はスキップして次の特徴量に移ります。

（結構長いです。ここはあくまで計算量を減らすテクニックなので次の章に移ってもいいです。）

while (f_i > n_total_constants and
      # Stop early if remaining features
      # are constant
      (n_visited_features < max_features or
      # At least one drawn features must be non constant                 
       n_visited_features <= n_found_constants + 
                             n_drawn_constants)):            
       
       n_visited_features +=1
       
       f_j = rand_int(n_drawn_constants, f_i - n_found_constants,                           
                      random_state)
       if f_j < n_known_constants:
       # f_j in the interval [n_drawn_constants, 
       # n_known_constants[
           features[n_drawn_constants], features[f_j] = 
           features[f_j], features[n_drawn_constants]                
           n_drawn_constants +=1　//14
           
       else:
           # f_j in the interval [n_known_constants, f_i - 
           # n_found_constants[
           f_j += n_found_constants
           # f_j in the interval [n_total_constants, f_i[                
           current.feature = features[f_j]
           # Sort samples along that feature; by
           # copying the values into an array and
           # sorting the array in a manner which utilizes the 
           cache more effectively.
           
           for i in range(start, end):    
               Xf[i] =self.X[samples[i], current.feature]                
           
           //15 データを取り出し小さい順にソート    
           sort(Xf + start, samples + start, end - start)
          
           if Xf[end -1] <= Xf[start] + FEATURE_THRESHOLD:                     
               features[f_j], features[n_total_constants] = 
               features[n_total_constants], features[f_j]                    
               // 16 Xが定数
               
               n_found_constants +=1                    
               n_total_constants +=1
               
           else:
               f_i -=1
                 
               // 17 Xが定数でない
               features[f_i], features[f_j] = 
               features[f_j], features[f_i]

探索中、特徴量は以下のように並んでいます。

: n_drawn_constants ・・・ 前回の調査で定数だと分かっている、かつ今回すでに当たったもの

n_drwan_constants: n_known_constants ・・・ 前回の調査ですでに定数だと分かっているもの

n_known_constants: n_found_constants ・・・ 今回の調査で定数だと分かったもの

n_found_constants: fi ・・・前回の調査では定数ではなかったもので、今回まだ当たっていないもの

fi: ・・・ 今回すでに当り、定数でないと分かったもの

各特徴量を当たったら、この規則が成り立つように並び替えます。

まず（省略してしまいましたが）fiは最後の特徴量からスタートします。n_drawn_constantsからfi-n_found_constantsの範囲で（未調査の範囲）特徴量fjを選びます。

・fjがn_known_constantsより小さい場合は、その特徴量は定数なのでfjをn_drawn_constantsの次に移動して次のループへ行きます。（14）

fj<n_known_constantsの場合

・それ以外の場合はデータXのfjの列を取り出し、値の小さい順に並べます。(15)

　⭐︎ Xの最大が最小とほぼ同じならば定数とみなしてn_found_constantsに追加し、次のル

　　ープへ行きます。(16)

　⭐︎ Xが定数でなければfiを1つ前に進めてfjと入れ替えます。その後、分割位置の探索に入

　　ります。(17)

n_known_constants<fjの場合

分割位置のループ

self.criterion.reset() 
p = start

while p < end:
    while (p +1< end and
           Xf[p +1] <= Xf[p] + FEATURE_THRESHOLD):                              
        //18 Xの値が変わらない場合はスキップ
        p +=1
                                   
    p +=1
    
    if p < end:
        current.pos = p
        
        # Reject if min_samples_leaf is not guaranteed
        if (((current.pos - start) < min_samples_leaf) or                                    
            ((end - current.pos) < min_samples_leaf)):
            continue
            
        self.criterion.update(current.pos) //19 ターゲットの分割
        
        # Reject if min_weight_leaf is not satisfied
        if ((self.criterion.weighted_n_left < min_weight_leaf) or                                    
            (self.criterion.weighted_n_right < min_weight_leaf)):
            continue                            
            
        current_proxy_improvement 
        =self.criterion.proxy_impurity_improvement()
        //20 不純度の増減分を計算
        
        if current_proxy_improvement > best_proxy_improvement:                                
            best_proxy_improvement = current_proxy_improvement
            # sum of halves is used to avoid infinite value                                
            
            current.threshold = Xf[p -1] /2.0+ Xf[p] /2.0
            
            if ((current.threshold == Xf[p]) or                                    
                (current.threshold == INFINITY) or                                    
                (current.threshold ==-INFINITY)):                                    
                current.threshold = Xf[p -1]                                
                
            best = current  # copy　//21 最適な分割

やっと分割位置（ここでは変数p）の探索が始まります。

まず、前後でXの値がほとんど変わらない場合は分割位置とはならないので飛ばします（18）。

19で選んだ分割位置でターゲット（y）を分割します。（詳細は後述）

20ではその分割によって不純度がどれだけ変わるかを計算しています。これまでの最適な分割位置よりも変化量が大きい場合はその分割位置を最適な分割位置として更新します（21）

if best.pos < end: 
    partition_end = end
    p = start
    
    //22 最適な分割位置で分割
    while p < partition_end:
        if self.X[samples[p], best.feature] <= best.threshold:                    
            p +=1
        
        else:
            partition_end -=1
            
            samples[p], samples[partition_end] = 
            samples[partition_end], samples[p]

     self.criterion.reset()
     self.criterion.update(best.pos)                        
     self.criterion.children_impurity(&best.impurity_left,
                                      &best.impurity_right)
            
     // 23 木全体で不純度を計算
     best.improvement =self.criterion.impurity_improvement(                
     impurity, best.impurity_left, best.impurity_right)
     
     # Respect invariant for constant features: the original 
     order of
     # element in features[:n_known_constants] must be preserved 
     for sibling
     # and child nodes
     memcpy(features, constant_features, sizeof(SIZE_t) * 
            n_known_constants)
            
     # Copy newly found constant features           
     memcpy(constant_features + n_known_constants,               
            features + n_known_constants,
             sizeof(SIZE_t) * n_found_constants)
            
     # Return values
     split[0] = best
     n_constant_features[0] = n_total_constants
     return0

ループが終わったら不純度の変化量が最大になる分割位置でデータを分割します（22）。

そして、その分割位置で木全体の不純度を計算します（23）。

これでデータの分割処理は終了です。以下ではターゲットを分割し、不純度を計算する過程を見ていきます。

ターゲットの分割

19のself.criterion.updateの中を見ていきます。

criterionはCriterion型のインスタンスです。Criterionクラスからは複数のクラスが派生していますが、今回はEntropyクラスを例にしたいと思います。（EntropyはClassificationCriterionクラスから派生します）

sklearn/tree/_criterion.pyx

cdef class ClassificationCriterion(Criterion):
"""Abstract criterion for classification."""

cdef class Entropy(ClassificationCriterion):
r"""Cross Entropy impurity criterion.    

This handles cases where the target is a classification taking values    0, 1, ... K-2, K-1. If node m represents a region Rm with Nm observations,    then let

    count_k = 1 / Nm \sum_{x_i in Rm} I(yi = k)
    
be the proportion of class k observations in node m.

The cross-entropy is then defined as

    cross-entropy = -\sum_{k=0}^{K-1} count_k log(count_k)
"""

updateメソッドはClassificationCriterionクラスのものをそのまま使います。

cdef int update(self, SIZE_t new_pos) nogil except-1:

    //略　変数の定義など

    if (new_pos - pos) <= (end - new_pos):
        for p in range(pos, new_pos):
            i = samples[p]
        
            if sample_weight !=NULL:
                w = sample_weight[i]
                
            for k in range(self.n_outputs):                     
                label_index = k *self.sum_stride +
                              <SIZE_t>self.y[i, k]                    
                sum_left[label_index] += w
            
            self.weighted_n_left += w
            
    else:
        self.reverse_reset()

        for p in range(end -1, new_pos -1, -1):
            i = samples[p]
            
            if sample_weight !=NULL:
                w = sample_weight[i]
                
            for k inrange(self.n_outputs):                    
                label_index = k *self.sum_stride +
                              <SIZE_t>self.y[i, k]                    
                sum_left[label_index] -= w
                
            self.weighted_n_left -= w
            
    # Update right part statistics
    self.weighted_n_right =self.weighted_n_node_samples -
                           self.weighted_n_left
    for k in range(self.n_outputs):
        for c in range(n_classes[k]):
            sum_right[c] = sum_total[c] - sum_left[c]            
            
        sum_right +=self.sum_stride
        sum_left +=self.sum_stride
        sum_total +=self.sum_stride
        
    self.pos = new_pos
    return0

n_outputsはyの次元数（つまりターゲットデータの種類数）です。ここではn_output=1として考えます。（多次元の場合も同じループが回るだけです）

sum_stridesはカテゴリ数です（多次元の場合はその最大値）。

sum_leftはカテゴリ数のサイズを持つ配列で、各要素は分割より左（つまりインデックスが小さい）のデータのうち、yがその値になる個数を表しています。sum_rightは右側、sum_totalはデータ全体について同様です。

weighted_n_left、weighted_n_right、weighted_n_node_samplesは左側、右側、全体のデータ数です。

ここでは一から計算するのではなくて、前回からの差分を計算しています。if文内は左が増えた場合、else文内は右が増えた場合です。

不純度の増減の計算

ターゲットを分割したのでそれに基づき不純度の増減を計算します。Entropyクラスのchildren_impurityメソッドです。

cdef void children_impurity(self, double* impurity_left,double* impurity_right) nogil:

//略　変数の定義など

    for k in range(self.n_outputs):
        for c in range(n_classes[k]):
            count_k = sum_left[c]
            if count_k >0.0:
                count_k /=self.weighted_n_left                     
                entropy_left -= count_k * log(count_k)  
                
            count_k = sum_right[c]
            if count_k >0.0:
                count_k /=self.weighted_n_right                    
                entropy_right -= count_k * log(count_k)            
                
        sum_left +=self.sum_stride 
        sum_right +=self.sum_stride
        
    impurity_left[0] = entropy_left /self.n_outputs        
    impurity_right[0] = entropy_right /self.n_outputs

これは定義通りエントロピーを計算しています。

参考文献

https://github.com/scikit-learn/scikit-learn/tree/master/sklearn/tree

決定木(Decision Tree)を理解して文書分類を行う

[入門]初心者の初心者による初心者のための決定木分析

最後に

機械学習の本質に関わる部分なのか計算量を減らすためのテクニックなのか判別するのが若干大変でした。

コメントとか変数の定義とか一部の処理を省略してしまったので気になる方は是非ご自身で読んでみてください。

あと、pythonだと思っているとポインタが出てきたり、C++だと思っていると定義せずに変数が使われたりとなかなか難しかったです。