C++實現GBDT算法及優化歷程

時間 2019-11-22

標籤 c++ 實現 gbdt 算法優化歷程欄目 C&C++ 简体版

原文原文鏈接

實際上是一個課程做業，要求實現 GBDT 算法。在實現的過程當中參考了不少資料，也作了不少優化，以爲收穫很大，所以把開發的過程也記錄了下來。node

源代碼在 GitHub。c++

1. 構建與使用

1.1 構建

Windows: 使用 Visual Studio 2017 打開解決方案並生成便可。
Linux: 根目錄提供了 makefile 文件，使用 make 編譯便可，須要 gcc >= 5.4.0

1.2 使用

用法：boost <config_file> <train_file> <test_file> <predict_dest>git

接受 LibSVM 格式的訓練數據輸入，以下每行表明一個訓練樣本：github

<label> <feature-index>:<feature-value> <feature-index>:<feature-value> <feature-index>:<feature-value>
複製代碼

用於預測的數據輸入和訓練數據相似：算法

<id> <feature-index>:<feature-value> <feature-index>:<feature-value> <feature-index>:<feature-value>
複製代碼

目前只支持二分類問題安全

<config_file> 指定訓練參數：性能優化

eta = 1.                 # shrinkage rate
gamma = 0.               # minimum gain required to split a node
maxDepth = 6             # max depth allowed
minChildWeight = 1       # minimum allowed size for a node to be splitted
rounds = 1               # REQUIRED. number of subtrees
subsample = 1.           # subsampling ratio for each tree
colsampleByTree = 1.     # tree-wise feature subsampling ratio
maxThreads = 1;          # max running threads
features;                # REQUIRED. number of features
validateSize = .2        # if greater than 0, input data will be split into two sets and used for training and validation repectively
複製代碼

2. 算法原理

GBDT 的核心能夠分紅兩部分，分別是 Gradient Boosting 和 Decision Tree:bash

Decision Tree : GBDT 的基分類器，經過劃分輸入樣本的特徵使得落在相同特徵的樣本擁有大體相同的 label。因爲在 GBDT 中須要對若干不一樣的 Decision Tree 的結果進行綜合，所以通常採用的是 Regression Tree （迴歸樹）而不是 Classification Tree （分類樹）。
Gradient Boosting: 迭代式的集成算法，每一棵決策樹的學習目標 y 都是以前全部樹的結論和的殘差（即梯度方向），也即 $y_i=y- \sum_{j=0}^{i-1}{\hat{y_j}}$ 。

3. 實現與優化歷程

各個部分的實現均通過若干次「第一版實現 - 性能 profiling - 優化獲得下一版代碼」的迭代。其中，性能 profiling 部分，使用的是 Visual Studio 2017 的「性能探查器」功能，在進行性能 profile 以前均使用 release 模式編譯（打開/O2 /Oi優化選項）。多線程

3.1 數據處理

選擇的輸入文件數據格式是 Libsvm 的格式，格式以下：dom

<label> <feature-index>:<feature-value> <feature-index>:<feature-value>
複製代碼

能夠看到這種格式自然適合用來表示稀疏的數據集，但在實現過程當中，爲了簡單起見以及 cache 性能，我經過將空值填充爲 0 轉化爲密集矩陣形式存儲。代價是內存佔用會相對高許多。

3.1.1 第一版

最初並無作什麼優化，採用的是以下的簡單流程：

文件按行讀取
對於每一行內容，先轉成 std::stringstream，再從中解析出相應的數據。

核心代碼以下：

ifstream in(path);
string line;
while (getline(in, line)) {
    auto item = parseLibSVMLine(move(line), featureCount); // { label, vector }
    x.push_back(move(item.first));
    y.push_back(item.second);
}

/* in parseLibSVMLine */
stringstream ss(line);
ss >> label;
while (ss) {
    char _;
    ss >> index >> _ >> value;
    values[index - 1] = value;
}
複製代碼

profile 結果：

能夠看到，主要的耗時在於將一行字符串解析成咱們須要的 label + vector 數據這一過程當中，進一步分析：

所以得知主要問題在於字符串解析部分。此時懷疑是 std::stringstream 的實現爲了線程安全、錯誤檢查等功能犧牲了性能，所以考慮使用 cstdio 中的實現。

3.1.2 改進

將 parseLibSVMLine 的實現重寫，使用cstdio 中的sscanf 代替了 std::stringstream：

int lastp = -1;
for (size_t p = 0; p < line.length(); p++) {
    if (isspace(line[p]) || p == line.length() - 1) {
        if (lastp == -1) {
            sscanf(line.c_str(), "%zu", &label);
        }
        else {
            sscanf(line.c_str() + lastp, "%zu:%lf", &index, &value);
            values[index - 1] = value;
        }
        lastp = int(p + 1);
    }
}
複製代碼

profile 結果：

能夠看到，雖然 parse 部分仍然是計算的熱點，但這部分的計算量顯著降低（53823 -> 23181），讀取完整個數據集的是時間減小了 50% 以上。

3.1.3 最終版

顯然，在數據集中，每一行之間的解析任務都是相互獨立的，所以能夠在一次性讀入整個文件並按行劃分數據後，對數據的解析進行並行化：

string content;
getline(ifstream(path), content, '\0');
stringstream in(move(content));

vector<string> lines;
string line;
while (getline(in, line)) lines.push_back(move(line));

#pragma omp parallel for
for (int i = 0; i < lines.size(); i++) {
    auto item = parseLibSVMLine(move(lines[i]), featureCount);
    #pragma omp critical
    {
        x.push_back(move(item.first));
        y.push_back(item.second);
    }
}
複製代碼

根據 profile 結果，進行並行化後，性能提高了約 25%。CPU 峯值佔用率從 15% 上升到了 70%。能夠發現性能的提高並無 CPU 佔用率的提高高，緣由根據推測有如下兩點：

讀取文件的IO時間，在測試時使用的是 672MB 的數據集，所以光是讀取所有內容就佔了 50% 以上的時間
多線程同步的代價

3.2 決策樹生成

決策樹生成的過程採用的是 depth-first 深度優先的方式，即不斷向下劃分子樹直到遇到下面的終止條件之一：

達到限定的最大深度
劃分收益小於閾值
該節點中的樣本數小於閾值

大體代碼以下：

auto p = new RegressionTree();

// calculate value for prediction
p->average = calculateAverageY();

if (x.size() > nodeThres) {
    // try to split
    auto ret = findSplitPoint(x, y, index);
    if (ret.gain > 0 && maxDepth > 1) { // check splitablity
        // split points
        // ...
        // ...

        p->left = createNode(x, y, leftIndex, maxDepth - 1);
        p->right = createNode(x, y, rightIndex, maxDepth - 1);
    }
}
複製代碼

3.2.1 計算劃分點

在哪一個特徵的哪一個值上作劃分是決策樹生成過程當中最核心（也是最耗時）的部分。

問題描述以下：

對於數據集，咱們要找到特徵以及該特徵上的劃分點，知足 MSE (mean-square-error 均方偏差) 最小：

\begin{equation}\begin{split}
E(A,q)&=E_{left}+E_{right}\\
&=\frac{N_1}{N_1+N_2} \sum_{(x_i,y_i)\in D_1(A,q)}(y_i-c_1)^2+\frac{N_2}{N_1+N_2} \sum_{(x_i,y_i)\in D_2(A,q)}(y_i-c_2)^2
\end{split}\end{equation}

其中：

$c_i=\frac{1}{N_i}\sum{(x_i,y_i)\in D_i(A,q)}{y_i}$ ，即劃分後的樣本 label 均值。
，爲劃分後的子數據集。

等價地，若是用表示劃分收益：

其中，爲劃分前的 MSE: $E_p=\sum_{(x_i,y_i)\in D}{(y_i-c)^2}$ ， $c=\frac{1}{|D|} \sum_{(x_i,y_i)\in D}y_i$ 。

尋找最佳劃分點等價於尋找收益最高的劃分方案：

(A,q)={\arg\max}_{A,q}G(A,q)={\arg\min}_{A,q}E(A,q)

3.2.1.1 基於排序的實現

分析：

\begin{equation}\begin{split}
E_{left}(A,q)
&= \sum_{(x_i,y_i)\in D_1(A,q)}(y_i-c_1)^2\\
&= \sum_{(x_i,y_i)\in D_1(A,q)}y_i^2+N_1c_1^2-2c_1^2 \sum_{(x_i,y_i)\in D_1(A,q)}y_i
\end{split}\end{equation}

\begin{equation}\begin{split}
E_{right}(A,q)
&=\sum_{(x_i,y_i)\in D_2(A,q)}(y_i-c_2)^2\\
&=\sum_{(x_i,y_i)\in D_2(A,q)}y_i^2+N_2c_2^2-2c_2^2 \sum_{(x_i,y_i)\in D_2(A,q)}y_i
\end{split}\end{equation}

顯然， $E_{left}$ 與 $E_{right}$ 都只與分割點左邊（右邊）的部分和有關，所以能夠先排序、再從小到大枚舉分割點計算出全部分割狀況的收益，對於每一個特徵，時間複雜度均爲 $O(n\log n)+O(n)=O(n)$ 。

代碼以下：

for (size_t featureIndex = 0; featureIndex < x.front().size(); featureIndex++) {
    vector<pair<size_t, double>> v(index.size());

    for (size_t i = 0; i < index.size(); i++) {
        auto ind = index[i];
        v[i].first = ind;
        v[i].second = x[ind][featureIndex];
    }

    // sorting
    tuple<size_t, double, double> tup;
    sort(v.begin(), v.end(), [](const auto &l, const auto &r) {
        return l.second < r.second;
    });

    // maintaining sums of y_i and y_i^2 in both left and right part
    double wholeErr, leftErr, rightErr;
    double wholeSum = 0, leftSum, rightSum;
    double wholePowSum = 0, leftPowSum, rightPowSum;
    for (const auto &t : v) {
        wholeSum += y[t.first];
        wholePowSum += pow(y[t.first], 2);
    }
    wholeErr = calculateError(index.size(), wholeSum, wholePowSum);

    leftSum = leftPowSum = 0;
    rightSum = wholeSum;
    rightPowSum = wholePowSum;
    for (size_t i = 0; i + 1 < index.size(); i++) {
        auto label = y[v[i].first];

        leftSum += label;
        rightSum -= label;
        leftPowSum += pow(label, 2);
        rightPowSum -= pow(label, 2);

        if (y[v[i].first] == y[v[i + 1].first]) continue; // same label with next, not splitable
        if (v[i].second == v[i + 1].second) continue; // same value, not splitable

        leftErr = calculateError(i + 1, leftSum, leftPowSum);
        rightErr = calculateError(index.size() - i - 1, rightSum, rightPowSum);

        // calculate error gain
        double gain = wholeErr - ((i + 1) * leftErr / index.size() + (index.size() - i - 1) * rightErr / index.size());
        if (gain > bestGain) {
            bestGain = gain;
            bestSplit = (v[i].second + v[i + 1].second) / 2;
            bestFeature = featureIndex;
        }
    }
}
複製代碼

profile 結果：

能夠看到， sorting 以及 sorting 以前的數據準備部分佔了很大一部分時間。

3.2.1.2 基於採樣分桶的實現

因爲以前基於排序的實現耗時較大，所以考慮換一種方法。後來翻 LightGBM 的優化方案，在參考文獻[^1]裏看到一個叫作 Sampling the Splitting points (SS) 的方法，比起 LightGBM 的方案， SS 方法更加容易實現。

SS 方法描述以下：

對於個亂序的數值，咱們先從中隨機採樣個樣本，將其排序後再等距採樣個樣本，以這個樣本做爲個桶的分割點。文獻中指出，若是，那麼有很高的機率能保證分到個桶中的樣本數量都接近 $\frac{n}{q}$ ，也就是接近等分。

採用這種方法，只須要的時間採樣出個桶、 $O(N\log q)$ 的時間來將全部樣本分配到不一樣的桶中。

在劃分桶以後，咱們只選擇桶的分割點做爲節點分割點的候選，所以只須要對代碼稍做改動便可在的時間內找到最佳的分割點。所以對於每一個特徵，尋找最佳分割點的時間複雜度爲 $O(N\log q)$ 。

使用這種方法，雖然由於只考慮了以分桶邊界的值進行分割的狀況，不必定能找到最佳的分割，但由於 Boosting 方法其本質即是將許多「次優」決策樹進行結合，所以 SS 方法形成的損失是能夠接受的。

[^1]: Ranka, Sanjay, and V. Singh. 「CLOUDS: A decision tree classifier for large datasets.」 Proceedings of the 4th Knowledge Discovery and Data Mining Conference. 1998.

代碼以下（簡單選擇 $s=\sqrt{N}$ ， $q=\sqrt{s}$ ）：

雖然這樣的取值事實上會使 $q=\sqrt[4]N$ ，時間複雜度 $O(N\log q)=O(N\log N^{\frac{1}{4}})=O(N\log N)$ ，但在測試數據中 $N\sim10^6$ ， $q~\sim32$ 已經足夠小。而若繼續增大，則能夠簡單將設爲一個不大於的常數，影響不大。

/* in findSplitPoint */
size_t nSample = size_t(pow(num, .5)), nBin = size_t(pow(num, .25));
auto dividers = sampleBinsDivider(x, nSample, nBin);
vector<double> binSums(nBin, .0), binPowSums(nBin, .0);
vector<size_t> binSizes(nBin, 0);
for (int i = 0; i < num; i++) {
    auto value = getFeatureValue(featureIndex, i);
    auto into = decideWhichBin(dividers, value);
    auto label = y[i];
    binSums[into] += label;
    binPowSums[into] += pow(label, 2);
    binSizes[into]++;
}
複製代碼

另外：由於數據集中數據的分佈是十分稀疏的，也即有大部分都是 0 值，所以在 decideWhichBin 中若是加入對小於第一個分割點的特判，將能帶來約 20% 的時間減小：

size_t RegressionTree::decideWhichBin(const std::vector<double>& divider, double value) {
    if (divider.empty() || value <= divider.front()) return 0;
    if (value > divider.back()) return divider.size();
    auto it = lower_bound(divider.cbegin(), divider.cend(), value);
    return it - divider.cbegin();
  }
複製代碼

根據在相同數據集、相同參數的測試結果，使用 SS 方法的每輪迭代時間減小約 80%。

3.2.1.3 加入並行

顯然在尋找最佳的劃分方案時，在不一樣的特徵上尋找最佳劃分點的任務是相互獨立的，所以能夠在特徵層面實現並行：

#pragma omp parallel for
for (int i = 0; i < featureIndexes.size(); i++) {
  /* sampling, bining... */

  // for each divider
    #pragma omp critical
    if (gain > bestGain) {
        bestGain = gain;
        bestSplit = divider;
        bestFeature = featureIndex;
    }
}
複製代碼

加入並行優化後，CPU峯值佔用率從 15% 提高到 70%，每輪迭代時間減小約 60%。

3.2.2 節點生成

以 depth-first 的順序進行生成，直到遇到終止條件爲止：

auto p = new RegressionTree();
// calculate value for prediction
p->average = average(y);
if (index.size() > max<size_t>(1, config.minChildWeight)) { // if this node is big enough
    // try to split
    auto ret = findSplitPoint(xx, y, index, featureIndexes);
    if (ret.gain > config.gamma && leftDepth > 1) { // check splitablity
        /* split points ... */
        // start splitting
        if (leftIndex.size() != 0 && rightIndex.size() != 0) {
            p->isLeaf = false;
            p->featureIndex = ret.featureIndex;
            p->featureValue = ret.splitPoint;
            // recursively build left and right subtrees
            p->left = createNode(leftX, leftY, config, leftDepth - 1);
            p->right = createNode(rightX, rightY, config, leftDepth - 1);
        }
    }
}
複製代碼

3.3 預測

對於輸入的每一個樣本，根據相應樹節點的劃分條件不斷向下劃分直到遇到葉子節點爲止，此時以葉子結點中的訓練樣本的平均 label 做爲預測值：

if (isLeaf) return average;
if (r[featureIndex] <= featureValue) return left->predict(r);
else return right->predict(r);
複製代碼

顯然，不一樣樣本之間的預測任務是相互獨立的，所以能夠對樣本之間的預測作並行：

Data::DataColumn result(x.size());
#pragma omp parallel for
for (int i = 0; i < x.size(); i++) {
    result[i] = predict(x[i]);
}
複製代碼

3.4 Boosting

Boosting 部分相對比較簡單，只須要在每次生成一棵新的決策樹後維護一下殘差便可：

while (roundsLeft--) {
    auto subtree = RegressionTree::fit(xx, residual, config);
    auto pred = subtree->predict(x);
    pred *= config.eta; // shrinkage rate
    residual -= pred;
}
複製代碼

4. 其餘優化

4.1 性能優化

4.1.1 sample性能

本來在採樣分割點的時候使用的是 C++17 標準中的 std::sample：

vector<double> samples(s);
vector<size_t> sampleIndex(s);
sample(index.begin(), index.end(), sampleIndex.begin(), s, mt19937{ random_device{}() });
for (size_t i = 0; i < s; i++) samples[i] = v[sampleIndex[i]];
複製代碼

但從 profiling 結果來看， std::sample 有很嚴重的效率問題：

對比使用普通隨機抽樣的狀況：

vector<double> samples(s);
std::random_device rd;
auto gen = std::default_random_engine(rd());
std::uniform_int_distribution<size_t> dis(0, index.size() - 1);
for (size_t i = 0; i < s; i++) samples[i] = v[index[dis(gen)]];
複製代碼

能夠看到，不使用 std::sample 的話，每輪耗時能減小一半以上。

4.1.2 節點分割

在劃分左右子樹的數據時，若是直接劃分 X, Y 數據的話會須要比較多的內存操做時間，所以這裏選擇的作法是：X, Y 固定不變，採用劃分索引的方式進行，經過索引來得到在屬於該節點的樣本下標：

for (size_t i = 0; i < index.size(); i++) {
    auto ind = index[i];
    if (xx[ret.featureIndex][ind] <= ret.splitPoint) {
        leftIndex.push_back(ind); // to the left
    }
    else {
        rightIndex.push_back(ind); // to the right
    }
}
複製代碼

4.2 並行化

並行化的實現依靠的是 OpenMP，經過形如 #pragma omp parallel 的編譯宏指令實現。

在實現中，有如下幾處使用了並行：

輸入數據處理，見 2.1.3
尋找最佳分裂點，見 2.2.1.3
預測，見 2.3

4.3 cache性能優化

4.3.1 X重整

對於 LibSVM 格式的輸入數據來講，一個很直覺的存儲方式是以 $N \times dims$ 的形狀存儲。但縱觀整個算法，在訓練過程當中對數據的訪問都是固定維的連續訪問（即對全部樣本的某一特徵的讀取），這樣不連續的內存訪問會形成 cache 性能的降低。所以在訓練以前，我把 $N\times dims$ 的數據重整成了以特徵優先的 $dims\times N$ 形狀，這樣在訓練過程當中就只須要對 x[featureIndex] 進行連續讀取，對 cache 更友好。

4.3.2 索引排序

在 3.1.2 中提到，爲了減小內存的操做而使用索引的形式來傳遞樣本劃分信息。但在後來發現形成了性能的降低，通過排查發現是由於加入了 subsample 功能即「對於每棵子樹只使用訓練樣本的一部分進行訓練」。爲了實現這一功能，在生成初始索引的時候：

// generate subsample
auto sampleSize = size_t(y.size() * config.subsample);
Index index(sampleSize);
std::uniform_int_distribution<size_t> dis(0, y.size() - 1);
for (size_t i = 0; i < index.size(); i++) index[i] = dis(gen); // sample with replacement
複製代碼

獲得的索引是無序的，這也形成了形如 x[featureIndex][index[i]] 的遍歷讀取是亂序的、cache 不友好的。因而經過對生成的索引進行排序從而解決：

// generate subsample
auto sampleSize = size_t(y.size() * config.subsample);
Index index(sampleSize);
std::uniform_int_distribution<size_t> dis(0, y.size() - 1);
for (size_t i = 0; i < index.size(); i++) index[i] = dis(gen); // sample with replacement
sort(index.begin(), index.end()); // for cache
複製代碼

4.3.3 連續 label 值

相比於龐大的 X 數據，Y 只有一列，所以不採起索引方式，直接劃分紅左右子樹的 yLeft與 yRight，進一步提高 cache 友好度：

// during splitting
vector<size_t> leftIndex, rightIndex;
Data::DataColumn leftY, rightY;
for (size_t i = 0; i < index.size(); i++) {
    auto ind = index[i];
    if (xx[ret.featureIndex][ind] <= ret.splitPoint) {
        leftIndex.push_back(ind); // to the left
        leftY.push_back(y[i]);    // split y
    }
    else {
        rightIndex.push_back(ind); // to the right
        rightY.push_back(y[i]);    // split y
    }
}
複製代碼

5. 測試對比

5.1 性能

主要與 xgboost 對比。

測試環境：

i7-5700HQ + 16GB

Ubuntu 16.04 (Windows Subsystem for Linux)

g++ -std=c++17 -O3 -fopenmp -m64

訓練數據：

train: $1719691\times 201$

max-depth: 20

subsample = 95

colsample-by-tree = .93

本算法：
- 數據讀取：耗時較大，約需 22s
- 訓練：每輪平均耗時 32s
xgboost：
- 數據讀取：約需 4s
- 訓練：每輪平均耗時 40s

5.2 內存佔用

本算法：因爲是把數據做爲密集矩陣存儲，所以對內存的消耗較大。以 672MB 的訓練數據爲例，由於數據較爲稀疏，讀入內存處理以後的佔用空間膨脹到了 5GB。
xgboost：一樣的 672MB 訓練數據，運行時約消耗 1.4GB。

5.3 預測準確性

因爲本算法使用了 SS 方法，所以相同輪數下的預測準確率應該低於 xgboost，簡單測試以下：

本算法：
xgboost：

簡單測試的意思是測試時並無對提供給本算法的訓練參數進行調優，使用的是以下配置：

rounds = 5

features = 201

eta = .3

maxThreads = 16

gamma = 1e-4

minChildWeight = 10

maxDepth = 20

validateSize = 0

subsample = 0.9500

colsampleByTree = 0.9287

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。