C++ реализует алгоритм GBDT и процесс оптимизации

На самом деле это задание курса, в котором предлагается реализовать алгоритм GBDT. В процессе реализации я обращался к большому количеству материалов и делал много оптимизаций.Чувствую, что много выиграл, поэтому так же записал процесс разработки.

исходный код вGitHub.

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>
Примите ввод обучающих данных в формате LibSVM, где каждая следующая строка представляет обучающую выборку:
```
<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 можно разделить на две части: повышение градиента и дерево решений:

Дерево решений: базовый классификатор GBDT, который разделяет признаки входных выборок таким образом, чтобы выборки, относящиеся к одним и тем же характеристикам, имели примерно одинаковую метку. Поскольку в GBDT необходимо синтезировать результаты нескольких разных деревьев решений, дерево регрессии обычно используется вместо дерева классификации.
Повышение градиента: алгоритм итеративного ансамбля, цель обучения y каждого дерева решений - это остаток (т.е. направление градиента) итоговой суммы всех предыдущих деревьев, то есть $y_i=y- \sum_{j=0}^{i-1}{\hat{y_j}}$ .

3. Процесс внедрения и оптимизации

Реализация каждой части претерпела несколько итераций «реализация первой версии — профилирование производительности — оптимизация для получения следующей версии кода». Среди них часть профилирования производительности использует функцию «профилировщик производительности» Visual Studio 2017 и использует режим выпуска для компиляции перед профилем производительности (открыть/O2 /Oiварианты оптимизации).

3.1 Обработка данных

Выбранный формат данных входного файла соответствует формату Libsvm, и этот формат выглядит следующим образом:

<label> <feature-index>:<feature-value> <feature-index>:<feature-value>

Видно, что этот формат естественно подходит для представления разреженных наборов данных, но в процессе реализации, для простоты и производительности кеша, я преобразовал его в плотную матричную форму, заполнив нулевые значения 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;
}

результат профиля:

Видно, что основные временные затраты приходятся на процесс разбора строки строк на нужные нам метки + векторные данные, и дальнейший анализ:

Итак, я знаю, что основная проблема заключается в части синтаксического анализа строк. Сомнение в это время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);
    }
}

результат профиля:

Видно, что хотя часть синтаксического анализа по-прежнему является горячей точкой вычислений, количество вычислений этой части значительно сократилось (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);
    }
}

По результатам профиля, после распараллеливания производительность улучшается примерно на 25%. Пиковая загрузка ЦП увеличилась с 15% до 70%. Можно обнаружить, что улучшение производительности не так высоко, как уровень загрузки процессора.Предположим, что причины заключаются в следующих двух пунктах:

Время ввода-вывода при чтении файла, в тесте использовался набор данных 672 МБ, поэтому простое чтение всего содержимого составляло более 50% времени.
Стоимость многопоточной синхронизации

3.2 Генерация дерева решений

Процесс генерации дерева решений использует подход «сначала в глубину — сначала в глубину», то есть поддерево непрерывно делится вниз до тех пор, пока не встретится одно из следующих условий завершения:

Достичь заданной максимальной глубины
Выгода от деления меньше порога
Количество выборок в этом узле меньше порога

Примерный код выглядит следующим образом:

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 (среднеквадратичной ошибке):

\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}$ , то есть среднее значение меток разделенной выборки.
,представляет собой разделенный набор подданных.

Эквивалентно, если использоватьУказывает разделенный доход:

в, E_p Для 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;
        }
    }
}

результат профиля:

Видно, что сортировка и подготовка данных перед сортировкой занимают большую часть времени.

3.2.1.2 Реализация на основе корзины для выборки

Поскольку предыдущая реализация на основе сортировки требует много времени, рассматривается другой метод. позже превратилсяLightGBMСхема оптимизации см. в ссылке [^1] метод, называемыйSampling the Splitting points (SS)По сравнению со схемой LightGBM метод SS проще в реализации.

Метод СС описывается следующим образом:

зазначения не по порядку, мы сначала случайным образом выбираем из нихвыборки, сортировать их, а затем пробовать через равные промежутки q-1 образцы, с этим q-1 образцы какточка разделения для каждого ведра. В литературе отмечается, что если q << s , то с большой вероятностьюКоличество образцов в каждом ведре близко к $\frac{n}{q}$ , что почти равно.

Таким образом, только O(N) выбранное времябочки, $O(N\log q)$ время распределить все образцы по разным ведрам.

После разделения сегментов мы выбираем только точки разделения сегментов в качестве кандидатов на точки разделения узлов, поэтому нам нужно лишь немного изменить код, чтобы O(q) найти наилучшую точку разделения во времени. Таким образом, для каждой функции временная сложность поиска лучшей точки разделения равна $O(N\log q)$ .

Используя этот метод, хотя и рассматривается только случай разбиения по значению границы корзины, возможно, не удастся найти наилучшее разбиение, а потому, что суть метода Boosting заключается в объединении множества «субоптимальных» деревьев решений , поэтому потери, вызванные методом SS, приемлемы.

[^1]: Ранка, Санджай и В. Сингх, «ОБЛАКА: классификатор дерева решений для больших наборов данных», Материалы 4-й конференции по исследованию знаний и интеллектуальному анализу данных, 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();
  }

По результатам тестирования на том же наборе данных и тех же параметрах время каждой итерации методом СС сокращается примерно на 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;
    }
}

После добавления параллельной оптимизации пиковая загрузка ЦП увеличилась с 15% до 70%, а время каждой итерации сократилось примерно на 60%.

3.2.2 Генерация узлов

Строить в порядке глубины до тех пор, пока не будет выполнено условие завершения:

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 Прогноз

Для каждой входной выборки она непрерывно делится в соответствии с условиями деления соответствующего узла дерева, пока не встретится с конечным узлом, В это время в качестве прогнозируемого значения используется средняя метка обучающих выборок в листовом узле:

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

Часть повышения относительно проста, просто нужно поддерживать остаток после каждого создания нового дерева решений:

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 Пример исполнения

Первоначально при выборке точек разделения использовался стандарт 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]];

Но по результатам профилирования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 Оптимизация производительности кэша

4.3.1 X реформация

Для входных данных в формате LibSVM очень интуитивно понятным способом их хранения является $N \times dims$ хранение формы. Но глядя на весь алгоритм, доступ к данным в процессе обучения фиксируется. dims Непрерывный доступ к измерению (то есть считывание признака всех выборок), такой прерывистый доступ к памяти приведет к снижению производительности кэша. Поэтому перед тренировкой я положил $N\times dims$ Данные преобразуются в функции в первую очередь $dims\times N$ форме, так что во время тренировки толькоx[featureIndex]Непрерывное чтение более удобно для кэша.

4.3.2 Сортировка по индексу

Как упоминалось в 3.1.2, для сокращения операций с памятью для передачи информации о делении выборки используется форма индекса. Однако позже выяснилось, что производительность ухудшилась.После расследования выяснилось, что это произошло из-за того, что была добавлена функция подвыборки, то есть «для каждого поддерева для обучения используется только часть обучающей выборки». Чтобы выполнить эту функцию, при создании начального индекса:

// 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]]обходные чтения не по порядку и не подходят для кэширования. Итак, это решается сортировкой сгенерированного индекса:

// 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 Непрерывные значения меток

По сравнению с огромными данными X, Y имеет только один столбец, поэтому он не использует метод индекса и напрямую делится на левое и правое поддеревья.yLeftиyRight, чтобы еще больше повысить удобство кэширования:

// 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

Этот алгоритм:
- Чтение данных: занимает много времени, около 22 с
- Обучение: среднее время за раунд 32 с.
xgboost:
- Чтение данных: около 4 с
- Обучение: Среднее время за раунд 40 секунд.

5.2 Использование памяти

Этот алгоритм: Поскольку данные хранятся в виде плотной матрицы, потребление памяти велико. Возьмем в качестве примера 672 МБ обучающих данных, поскольку данные относительно разрежены, занятое пространство после считывания в память для обработки увеличилось до 5 ГБ.
xgboost: те же обучающие данные объемом 672 МБ потребляют около 1,4 ГБ во время выполнения.

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