source: trunk/sources/HeuristicLab.Algorithms.DataAnalysis/3.4/GradientBoostedTrees/RegressionTreeBuilder.cs @ 12700

Last change on this file since 12700 was 12700, checked in by gkronber, 5 years ago

#2261: copied GBT implementation from branch to trunk

File size: 20.7 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22
23using System;
24using System.Collections;
25using System.Collections.Generic;
26using System.Diagnostics;
27using System.Linq;
28using HeuristicLab.Core;
29using HeuristicLab.Problems.DataAnalysis;
30
31namespace HeuristicLab.Algorithms.DataAnalysis {
32  // This class implements a greedy decision tree learner which selects splits with the maximum reduction in sum of squared errors.
33  // The tree builder also tracks variable relevance metrics based on the splits and improvement after the split.
34  // The implementation is tuned for gradient boosting where multiple trees have to be calculated for the same training data
35  // each time with a different target vector. Vectors of idx to allow iteration of intput variables in sorted order are
36  // pre-calculated so that optimal thresholds for splits can be calculated in O(n) for each input variable.
37  // After each split the row idx are partitioned in a left an right part.
38  internal class RegressionTreeBuilder {
39    private readonly IRandom random;
40    private readonly IRegressionProblemData problemData;
41
42    private readonly int nCols;
43    private readonly double[][] x; // all training data (original order from problemData), x is constant
44    private double[] originalY; // the original target labels (from problemData), originalY is constant
45    private double[] curPred; // current predictions for originalY (in case we are using gradient boosting, otherwise = zeros), only necessary for line search
46
47    private double[] y; // training labels (original order from problemData), y can be changed
48
49    private Dictionary<string, double> sumImprovements; // for variable relevance calculation
50
51    private readonly string[] allowedVariables; // all variables in shuffled order
52    private Dictionary<string, int> varName2Index; // maps the variable names to column indexes
53    private int effectiveVars; // number of variables that are used from allowedVariables
54
55    private int effectiveRows; // number of rows that are used from
56    private readonly int[][] sortedIdxAll;
57    private readonly int[][] sortedIdx; // random selection from sortedIdxAll (for r < 1.0)
58
59    // helper arrays which are allocated to maximal necessary size only once in the ctor
60    private readonly int[] internalIdx, which, leftTmp, rightTmp;
61    private readonly double[] outx;
62    private readonly int[] outSortedIdx;
63
64    private RegressionTreeModel.TreeNode[] tree; // tree is represented as a flat array of nodes
65    private int curTreeNodeIdx; // the index where the next tree node is stored
66
67    // This class represents information about potential splits.
68    // For each node generated the best splitting variable and threshold as well as
69    // the improvement from the split are stored in a priority queue
70    private class PartitionSplits {
71      public int ParentNodeIdx { get; set; } // the idx of the leaf node representing this partition
72      public int StartIdx { get; set; } // the start idx of the partition
73      public int EndIndex { get; set; } // the end idx of the partition
74      public string SplittingVariable { get; set; } // the best splitting variable
75      public double SplittingThreshold { get; set; } // the best threshold
76      public double SplittingImprovement { get; set; } // the improvement of the split (for priority queue)
77    }
78
79    // this list hold partitions with the information about the best split (organized as a sorted queue)
80    private readonly IList<PartitionSplits> queue;
81
82    // prepare and allocate buffer variables in ctor
83    public RegressionTreeBuilder(IRegressionProblemData problemData, IRandom random) {
84      this.problemData = problemData;
85      this.random = random;
86
87      var rows = problemData.TrainingIndices.Count();
88
89      this.nCols = problemData.AllowedInputVariables.Count();
90
91      allowedVariables = problemData.AllowedInputVariables.ToArray();
92      varName2Index = new Dictionary<string, int>(allowedVariables.Length);
93      for (int i = 0; i < allowedVariables.Length; i++) varName2Index.Add(allowedVariables[i], i);
94
95      sortedIdxAll = new int[nCols][];
96      sortedIdx = new int[nCols][];
97      sumImprovements = new Dictionary<string, double>();
98      internalIdx = new int[rows];
99      which = new int[rows];
100      leftTmp = new int[rows];
101      rightTmp = new int[rows];
102      outx = new double[rows];
103      outSortedIdx = new int[rows];
104      queue = new List<PartitionSplits>(100);
105
106      x = new double[nCols][];
107      originalY = problemData.Dataset.GetDoubleValues(problemData.TargetVariable, problemData.TrainingIndices).ToArray();
108      y = new double[originalY.Length];
109      Array.Copy(originalY, y, y.Length); // copy values (originalY is fixed, y is changed in gradient boosting)
110      curPred = Enumerable.Repeat(0.0, y.Length).ToArray(); // zeros
111
112      int col = 0;
113      foreach (var inputVariable in problemData.AllowedInputVariables) {
114        x[col] = problemData.Dataset.GetDoubleValues(inputVariable, problemData.TrainingIndices).ToArray();
115        sortedIdxAll[col] = Enumerable.Range(0, rows).OrderBy(r => x[col][r]).ToArray();
116        sortedIdx[col] = new int[rows];
117        col++;
118      }
119    }
120
121    // simple API produces a single regression tree optimizing sum of squared errors
122    // this can be used if only a simple regression tree should be produced
123    // for a set of trees use the method CreateRegressionTreeForGradientBoosting below
124    //
125    // r and m work in the same way as for alglib random forest
126    // r is fraction of rows to use for training
127    // m is fraction of variables to use for training
128    public IRegressionModel CreateRegressionTree(int maxSize, double r = 0.5, double m = 0.5) {
129      // subtract mean of y first
130      var yAvg = y.Average();
131      for (int i = 0; i < y.Length; i++) y[i] -= yAvg;
132
133      var seLoss = new SquaredErrorLoss();
134
135      var model = CreateRegressionTreeForGradientBoosting(y, curPred, maxSize, problemData.TrainingIndices.ToArray(), seLoss, r, m);
136
137      return new GradientBoostedTreesModel(new[] { new ConstantRegressionModel(yAvg), model }, new[] { 1.0, 1.0 });
138    }
139
140    // specific interface that allows to specify the target labels and the training rows which is necessary when for gradient boosted trees
141    public IRegressionModel CreateRegressionTreeForGradientBoosting(double[] y, double[] curPred, int maxSize, int[] idx, ILossFunction lossFunction, double r = 0.5, double m = 0.5) {
142      Debug.Assert(maxSize > 0);
143      Debug.Assert(r > 0);
144      Debug.Assert(r <= 1.0);
145      Debug.Assert(y.Count() == this.y.Length);
146      Debug.Assert(m > 0);
147      Debug.Assert(m <= 1.0);
148
149      // y and curPred are changed in gradient boosting
150      this.y = y;
151      this.curPred = curPred;
152
153      // shuffle row idx
154      HeuristicLab.Random.ListExtensions.ShuffleInPlace(idx, random);
155
156      int nRows = idx.Count();
157
158      // shuffle variable idx
159      HeuristicLab.Random.ListExtensions.ShuffleInPlace(allowedVariables, random);
160
161      // only select a part of the rows and columns randomly
162      effectiveRows = (int)Math.Ceiling(nRows * r);
163      effectiveVars = (int)Math.Ceiling(nCols * m);
164
165      // the which array is used for partitioing row idxs 
166      Array.Clear(which, 0, which.Length);
167
168      // mark selected rows
169      for (int row = 0; row < effectiveRows; row++) {
170        which[idx[row]] = 1; // we use the which vector as a temporary variable here
171        internalIdx[row] = idx[row];
172      }
173
174      for (int col = 0; col < nCols; col++) {
175        int i = 0;
176        for (int row = 0; row < nRows; row++) {
177          if (which[sortedIdxAll[col][row]] > 0) {
178            Debug.Assert(i < effectiveRows);
179            sortedIdx[col][i] = sortedIdxAll[col][row];
180            i++;
181          }
182        }
183      }
184
185      this.tree = new RegressionTreeModel.TreeNode[maxSize];
186      this.queue.Clear();
187      this.curTreeNodeIdx = 0;
188
189      // start out with only one leaf node (constant prediction)
190      // and calculate the best split for this root node and enqueue it into a queue sorted by improvement throught the split
191      // start and end idx are inclusive
192      CreateLeafNode(0, effectiveRows - 1, lossFunction);
193
194      // process the priority queue to complete the tree
195      CreateRegressionTreeFromQueue(maxSize, lossFunction);
196
197      return new RegressionTreeModel(tree.ToArray());
198    }
199
200
201    // processes potential splits from the queue as long as splits are left and the maximum size of the tree is not reached
202    private void CreateRegressionTreeFromQueue(int maxNodes, ILossFunction lossFunction) {
203      while (queue.Any() && curTreeNodeIdx + 1 < maxNodes) { // two nodes are created in each loop
204        var f = queue[queue.Count - 1]; // last element has the largest improvement
205        queue.RemoveAt(queue.Count - 1);
206
207        var startIdx = f.StartIdx;
208        var endIdx = f.EndIndex;
209
210        Debug.Assert(endIdx - startIdx >= 0);
211        Debug.Assert(startIdx >= 0);
212        Debug.Assert(endIdx < internalIdx.Length);
213
214        // split partition into left and right
215        int splitIdx;
216        SplitPartition(f.StartIdx, f.EndIndex, f.SplittingVariable, f.SplittingThreshold, out splitIdx);
217        Debug.Assert(splitIdx + 1 <= endIdx);
218        Debug.Assert(startIdx <= splitIdx);
219
220        // create two leaf nodes (and enqueue best splits for both)
221        var leftTreeIdx = CreateLeafNode(startIdx, splitIdx, lossFunction);
222        var rightTreeIdx = CreateLeafNode(splitIdx + 1, endIdx, lossFunction);
223
224        // overwrite existing leaf node with an internal node
225        tree[f.ParentNodeIdx] = new RegressionTreeModel.TreeNode(f.SplittingVariable, f.SplittingThreshold, leftTreeIdx, rightTreeIdx);
226      }
227    }
228
229
230    // returns the index of the newly created tree node
231    private int CreateLeafNode(int startIdx, int endIdx, ILossFunction lossFunction) {
232      // write a leaf node
233      var val = lossFunction.LineSearch(originalY, curPred, internalIdx, startIdx, endIdx);
234      tree[curTreeNodeIdx] = new RegressionTreeModel.TreeNode(RegressionTreeModel.TreeNode.NO_VARIABLE, val);
235
236      EnqueuePartitionSplit(curTreeNodeIdx, startIdx, endIdx);
237      curTreeNodeIdx++;
238      return curTreeNodeIdx - 1;
239    }
240
241
242    // calculates the optimal split for the partition [startIdx .. endIdx] (inclusive)
243    // which is represented by the leaf node with the specified nodeIdx
244    private void EnqueuePartitionSplit(int nodeIdx, int startIdx, int endIdx) {
245      double threshold, improvement;
246      string bestVariableName;
247      // only enqueue a new split if there are at least 2 rows left and a split is possible
248      if (startIdx < endIdx &&
249        FindBestVariableAndThreshold(startIdx, endIdx, out threshold, out bestVariableName, out improvement)) {
250        var split = new PartitionSplits() {
251          ParentNodeIdx = nodeIdx,
252          StartIdx = startIdx,
253          EndIndex = endIdx,
254          SplittingThreshold = threshold,
255          SplittingVariable = bestVariableName
256        };
257        InsertSortedQueue(split);
258      }
259    }
260
261
262    // routine for splitting a partition of rows stored in internalIdx between startIdx and endIdx into
263    // a left partition and a right partition using the given splittingVariable and threshold
264    // the splitIdx is the last index of the left partition
265    // splitIdx + 1 is the first index of the right partition
266    // startIdx and endIdx are inclusive
267    private void SplitPartition(int startIdx, int endIdx, string splittingVar, double threshold, out int splitIdx) {
268      int bestVarIdx = varName2Index[splittingVar];
269      // split - two pass
270
271      // store which index goes into which partition
272      for (int k = startIdx; k <= endIdx; k++) {
273        if (x[bestVarIdx][internalIdx[k]] <= threshold)
274          which[internalIdx[k]] = -1; // left partition
275        else
276          which[internalIdx[k]] = 1; // right partition
277      }
278
279      // partition sortedIdx for each variable
280      int i;
281      int j;
282      for (int col = 0; col < nCols; col++) {
283        i = 0;
284        j = 0;
285        int k;
286        for (k = startIdx; k <= endIdx; k++) {
287          Debug.Assert(Math.Abs(which[sortedIdx[col][k]]) == 1);
288
289          if (which[sortedIdx[col][k]] < 0) {
290            leftTmp[i++] = sortedIdx[col][k];
291          } else {
292            rightTmp[j++] = sortedIdx[col][k];
293          }
294        }
295        Debug.Assert(i > 0); // at least on element in the left partition
296        Debug.Assert(j > 0); // at least one element in the right partition
297        Debug.Assert(i + j == endIdx - startIdx + 1);
298        k = startIdx;
299        for (int l = 0; l < i; l++) sortedIdx[col][k++] = leftTmp[l];
300        for (int l = 0; l < j; l++) sortedIdx[col][k++] = rightTmp[l];
301      }
302
303      // partition row indices
304      i = startIdx;
305      j = endIdx;
306      while (i <= j) {
307        Debug.Assert(Math.Abs(which[internalIdx[i]]) == 1);
308        Debug.Assert(Math.Abs(which[internalIdx[j]]) == 1);
309        if (which[internalIdx[i]] < 0) i++;
310        else if (which[internalIdx[j]] > 0) j--;
311        else {
312          Debug.Assert(which[internalIdx[i]] > 0);
313          Debug.Assert(which[internalIdx[j]] < 0);
314          // swap
315          int tmp = internalIdx[i];
316          internalIdx[i] = internalIdx[j];
317          internalIdx[j] = tmp;
318          i++;
319          j--;
320        }
321      }
322      Debug.Assert(j + 1 == i);
323      Debug.Assert(i <= endIdx);
324      Debug.Assert(startIdx <= j);
325
326      splitIdx = j;
327    }
328
329    private bool FindBestVariableAndThreshold(int startIdx, int endIdx, out double threshold, out string bestVar, out double improvement) {
330      Debug.Assert(startIdx < endIdx + 1); // at least 2 elements
331
332      int rows = endIdx - startIdx + 1;
333      Debug.Assert(rows >= 2);
334
335      double sumY = 0.0;
336      for (int i = startIdx; i <= endIdx; i++) {
337        sumY += y[internalIdx[i]];
338      }
339
340      // see description of calculation in FindBestThreshold
341      double bestImprovement = 1.0 / rows * sumY * sumY; // any improvement must be larger than this baseline
342      double bestThreshold = double.PositiveInfinity;
343      bestVar = RegressionTreeModel.TreeNode.NO_VARIABLE;
344
345      for (int col = 0; col < effectiveVars; col++) {
346        // sort values for variable to prepare for threshold selection
347        var curVariable = allowedVariables[col];
348        var curVariableIdx = varName2Index[curVariable];
349        for (int i = startIdx; i <= endIdx; i++) {
350          var sortedI = sortedIdx[curVariableIdx][i];
351          outSortedIdx[i - startIdx] = sortedI;
352          outx[i - startIdx] = x[curVariableIdx][sortedI];
353        }
354
355        double curImprovement;
356        double curThreshold;
357        FindBestThreshold(outx, outSortedIdx, rows, y, sumY, out curThreshold, out curImprovement);
358
359        if (curImprovement > bestImprovement) {
360          bestImprovement = curImprovement;
361          bestThreshold = curThreshold;
362          bestVar = allowedVariables[col];
363        }
364      }
365      if (bestVar == RegressionTreeModel.TreeNode.NO_VARIABLE) {
366        // not successfull
367        threshold = double.PositiveInfinity;
368        improvement = double.NegativeInfinity;
369        return false;
370      } else {
371        UpdateVariableRelevance(bestVar, sumY, bestImprovement, rows);
372        improvement = bestImprovement;
373        threshold = bestThreshold;
374        return true;
375      }
376    }
377
378    // x [0..N-1] contains rows sorted values in the range from [0..rows-1]
379    // sortedIdx [0..N-1] contains the idx of the values in x in the original dataset in the range from [0..rows-1]
380    // rows specifies the number of valid entries in x and sortedIdx
381    // y [0..N-1] contains the target values in original sorting order
382    // sumY is y.Sum()
383    //
384    // the routine returns the best threshold (x[i] + x[i+1]) / 2 for i = [0 .. rows-2] by calculating the reduction in squared error
385    // additionally the reduction in squared error is returned in bestImprovement
386    // if all elements of x are equal the routing fails to produce a threshold
387    private static void FindBestThreshold(double[] x, int[] sortedIdx, int rows, double[] y, double sumY, out double bestThreshold, out double bestImprovement) {
388      Debug.Assert(rows >= 2);
389
390      double sl = 0.0;
391      double sr = sumY;
392      double nl = 0.0;
393      double nr = rows;
394
395      bestImprovement = 1.0 / rows * sumY * sumY; // this is the baseline for the improvement
396      bestThreshold = double.NegativeInfinity;
397      // for all thresholds
398      // if we have n rows there are n-1 possible splits
399      for (int i = 0; i < rows - 1; i++) {
400        sl += y[sortedIdx[i]];
401        sr -= y[sortedIdx[i]];
402
403        nl++;
404        nr--;
405        Debug.Assert(nl > 0);
406        Debug.Assert(nr > 0);
407
408        if (x[i] < x[i + 1]) { // don't try to split when two elements are equal
409
410          // goal is to find the split with leading to minimal total variance of left and right parts
411          // without partitioning the variance is var(y) = E(y²) - E(y)² 
412          //    = 1/n * sum(y²) - (1/n * sum(y))²
413          //      -------------   ---------------
414          //         constant       baseline for improvement
415          //
416          // if we split into right and left part the overall variance is the weigthed combination nl/n * var(y_l) + nr/n * var(y_r) 
417          //    = nl/n * (1/nl * sum(y_l²) - (1/nl * sum(y_l))²) + nr/n * (1/nr * sum(y_r²) - (1/nr * sum(y_r))²)
418          //    = 1/n * sum(y_l²) - 1/nl * 1/n * sum(y_l)² + 1/n * sum(y_r²) - 1/nr * 1/n * sum(y_r)²
419          //    = 1/n * (sum(y_l²) + sum(y_r²)) - 1/n * (sum(y_l)² / nl + sum(y_r)² / nr)
420          //    = 1/n * sum(y²) - 1/n * (sum(y_l)² / nl + sum(y_r)² / nr)
421          //      -------------
422          //       not changed by split (and the same for total variance without partitioning)
423          //
424          //   therefore we need to find the maximum value (sum(y_l)² / nl + sum(y_r)² / nr) (ignoring the factor 1/n)
425          //   and this value must be larger than 1/n * sum(y)² to be an improvement over no split
426
427          double curQuality = sl * sl / nl + sr * sr / nr;
428
429          if (curQuality > bestImprovement) {
430            bestThreshold = (x[i] + x[i + 1]) / 2.0;
431            bestImprovement = curQuality;
432          }
433        }
434      }
435
436      // if all elements where the same then no split can be found
437    }
438
439
440    private void UpdateVariableRelevance(string bestVar, double sumY, double bestImprovement, int rows) {
441      if (string.IsNullOrEmpty(bestVar)) return;
442      // update variable relevance
443      double baseLine = 1.0 / rows * sumY * sumY; // if best improvement is equal to baseline then the split had no effect
444
445      double delta = (bestImprovement - baseLine);
446      double v;
447      if (!sumImprovements.TryGetValue(bestVar, out v)) {
448        sumImprovements[bestVar] = delta;
449      }
450      sumImprovements[bestVar] = v + delta;
451    }
452
453    public IEnumerable<KeyValuePair<string, double>> GetVariableRelevance() {
454      // values are scaled: the most important variable has relevance = 100
455      double scaling = 100 / sumImprovements.Max(t => t.Value);
456      return
457        sumImprovements
458        .Select(t => new KeyValuePair<string, double>(t.Key, t.Value * scaling))
459        .OrderByDescending(t => t.Value);
460    }
461
462
463    // insert a new parition split (find insertion point and start at first element of the queue)
464    // elements are removed from the queue at the last position
465    // O(n), splits could be organized as a heap to improve runtime (see alglib tsort)
466    private void InsertSortedQueue(PartitionSplits split) {
467      // find insertion position
468      int i = 0;
469      while (i < queue.Count && queue[i].SplittingImprovement < split.SplittingImprovement) { i++; }
470
471      queue.Insert(i, split);
472    }
473  }
474}
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
Note: See TracBrowser for help on using the repository browser.