James D. McCaffrey

I fed my decision tree regression system, implemented using C#, to several AI systems and asked the AI to analyze it for correctness. I was quite impressed that the AI found a few rare edge cases where my code could fail, and showed me how to check for those edge cases.

Additionally, the AI pointed out that even though my implementation was functionally correct, it used nested loops which gave complexity of O(N^2). This is fine for datasets of up to about 2,000 items, but past that, training would slow to a crawl.

I knew this, but I also knew that writing a performant version is extremely difficult. I decided to bite the bullet, and use AI to write a performant version of decision tree regression.

The effort was every bit as difficult as I expected, and took about 16 hours, even with AI’s tireless help.

For my demo, I used one of my standard datasets It looks like:

-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
. . .

The data is synthetic. The first five values on each line are the predictors. The last value on each line is the target to predict. There are 200 training items and 40 test items.

The key parts of the output of my revised demo are:

Begin decision tree regression (performant version)

Loading synthetic train (200) and test (40) data
Done

First three train X:
 -0.1660  0.4406 -0.9998 -0.3953 -0.7065
  0.0776 -0.1616  0.3704 -0.5911  0.7562
 -0.9452  0.3409 -0.1654  0.1174 -0.7192

First three train y:
  0.4840
  0.1568
  0.8054

Setting maxDepth = 3
Setting minSamples = 2
Setting minLeaf = 18
Using default numSplitCols = -1

Creating and training tree
Done

Tree:
ID 0   | sc   0 | sv  -0.2102 | L   1 | R   2 | py   0.3493 | leaf F | rc 200
ID 1   | sc   4 | sv   0.1431 | L   3 | R   4 | py   0.5345 | leaf F | rc 75
ID 2   | sc   0 | sv   0.3915 | L   5 | R   6 | py   0.2382 | leaf F | rc 125
ID 3   | sc   0 | sv  -0.6553 | L   7 | R   8 | py   0.6358 | leaf F | rc 41
ID 4   | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.4123 | leaf T | rc 34
ID 5   | sc   4 | sv  -0.2987 | L  11 | R  12 | py   0.3032 | leaf F | rc 64
ID 6   | sc   2 | sv   0.3777 | L  13 | R  14 | py   0.1701 | leaf F | rc 61
ID 7   | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.6952 | leaf T | rc 23
ID 8   | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.5598 | leaf T | rc 18
ID 11  | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.4101 | leaf T | rc 18
ID 12  | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.2613 | leaf T | rc 46
ID 13  | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.1882 | leaf T | rc 39
ID 14  | sc  -1 | sv   0.0000 | L  -1 | R  -1 | py   0.1381 | leaf T | rc 22

Rows assoc with node [11]:
0 7 14 24 69 87 88 119 121 123 133 136 138 141 162 186 191 195

Evaluating model
Accuracy train (within 0.10) = 0.3750
Accuracy test (within 0.10) = 0.4750

MSE train = 0.0048
MSE test = 0.0054

Predicting for trainX[0] =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065
Predicted y = 0.4101

IF
column 0  >   -0.2102 AND
column 0  <=   0.3915 AND
column 4  <=  -0.2987 AND
THEN node [11] predicted = 0.4101

End demo

Each line of the tree shows node ID, split column, split value (aka threshold), left child idx, right child idx, predicted value for the node, a Boolean to tell if the node is a leaf node or not, the the row count of rows associated with the node. The explain output has a trailing “AND” with no condition, but I was too lazy to strip it away.

The diagram below shows how the prediction was arrived at.

The accuracy is low, which is expected. Decision trees are almost never used by themselves. Instead they are usually part of a collection — bagging tree regression, random forest regression, AdaBoost regression, gradient boost regression.

Good fun.

I love models of all kinds — math models, machine learning models, etc., etc., and model trains. Here’s a beautiful mine train model in HOn30 scale.

Demo program. Very long, extremely complex. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. (My blog editor chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

// full List storage with indexes (no pointers)
// iteration-based construction (no stack, no recursion)
// Nodes hold associated row idxs (interpretability)
// but rows can be deleted after training (ensembles)
// highly performant version via AI help

namespace DecisionTreeRegression
{
  internal class DecisionTreeRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin decision tree" +
        " regression (performant version) ");

      // 1. load data
      Console.WriteLine("\nLoading synthetic train" +
        " (200) and test (40) data");
      string trainFile = "..\\..\\..\\Data\\" +
        "synthetic_train_200.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4 };
      int colY = 5;
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        new int[] { colY }, ',', "#"));

      string testFile = "..\\..\\..\\Data\\" +
        "synthetic_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        new int[] { colY }, ',', "#"));
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three train X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(trainX[i], 4, 8);

      Console.WriteLine("\nFirst three train y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(trainY[i].ToString("F4").
          PadLeft(8));

      // 2. create and train/build tree
      // small tree for blog, article
      int maxDepth = 3; //
      int minSamples = 2;  // to try split
      int minLeaf = 18;    // after split
      int numSplitCols = -1;  // use all columns
      bool saveRows = true;

      Console.WriteLine("\nSetting maxDepth = " +
        maxDepth);
      Console.WriteLine("Setting minSamples = " +
        minSamples);
      Console.WriteLine("Setting minLeaf = " +
        minLeaf);
      Console.WriteLine("Using default numSplitCols = -1 ");

      Console.WriteLine("\nCreating and training tree ");
      DecisionTreeRegressor dtr =
        new DecisionTreeRegressor(maxDepth, minSamples,
        minLeaf, numSplitCols, saveRows, seed: 0);
      dtr.Train(trainX, trainY);
      Console.WriteLine("Done ");

      Console.WriteLine("\nTree: ");
      dtr.Display();

      Console.WriteLine("\nRows assoc with node [11]: ");
      for (int i = 0; i "lt" dtr.tree[11].rows.Count; ++i)
      {
        if (i "gt" 0 && i % 20 == 0) Console.WriteLine("");
        Console.Write(dtr.tree[11].rows[i] + " ");
      }
      Console.WriteLine("");

      // 3. evaluate model
      Console.WriteLine("\nEvaluating model ");
      double accTrain = dtr.Accuracy(trainX, trainY, 0.10);
      Console.WriteLine("Accuracy train (within 0.10) = " +
        accTrain.ToString("F4"));
      double accTest = dtr.Accuracy(testX, testY, 0.10);
      Console.WriteLine("Accuracy test (within 0.10) = " +
        accTest.ToString("F4"));

      double mseTrain = dtr.MSE(trainX, trainY);
      Console.WriteLine("\nMSE train = " +
        mseTrain.ToString("F4"));
      double mseTest = dtr.MSE(testX, testY);
      Console.WriteLine("MSE test = " +
        mseTest.ToString("F4"));

      // 4. use model
      Console.WriteLine("\nPredicting for trainX[0] = ");
      double[] x = trainX[0];
      VecShow(x, 4, 9);
      double predY = dtr.Predict(x);
      Console.WriteLine("Predicted y = " +
        predY.ToString("F4"));

      dtr.Explain(x);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    // ------------------------------------------------------
    // helpers for Main():
    //   MatLoad(), MatToVec(), VecShow().
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result =
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static double[] MatToVec(double[][] mat)
    {
      int nRows = mat.Length;
      int nCols = mat[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = mat[i][j];
      return result;
    }

    static void VecShow(double[] vec, int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

  } // class Program

  // ========================================================

  public class DecisionTreeRegressor
  {
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"Node"gt" tree;
    public Random rnd;
    public bool saveRows;
    public double[][] trainX;
    public double[] trainY;

    // ............................................

    public class Node
    {
      public int id;
      public int colIdx;
      public double thresh;
      public int left;
      public int right;
      public double value;
      public bool isLeaf;
      public List"lt"int"gt" rows;

      public Node()
      {
        this.id = -1;
        this.colIdx = -1;
        this.thresh = 0.0;
        this.left = -1;
        this.right = -1;
        this.value = 0.0;
        this.isLeaf = false;
        this.rows = null;
      }
    }

    // ............................................

    public DecisionTreeRegressor(int maxDepth = 3,
      int minSamples = 2, int minLeaf = 1,
      int numSplitCols = -1, bool saveRows = false,
      int seed = 0)
    {
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.saveRows = saveRows;
      this.tree = new List"lt"Node"gt"();

      int numNodes = (int)Math.Pow(2, (maxDepth + 1)) - 1;
      for (int i = 0; i "lt" numNodes; ++i)
        this.tree.Add(null);
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------
    // public: ctor(), Train(), Predict(), Explain(),
    //  Display(), Accuracy(), MSE()
    // private: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

    public void Train(double[][] trainX, double[] trainY)
    {
      this.trainX = trainX;
      this.trainY = trainY;

      // boundary IDs based on max allowed depth
      int maxID = (int)Math.Pow(2, (this.maxDepth + 1)) - 2;
      int maxStartID = (int)Math.Pow(2, this.maxDepth) - 1;

      // initialize root data rows
      List"lt"int"gt" allRows = new List"lt"int"gt"(this.trainX.Length);
      for (int i = 0; i "lt" this.trainX.Length; ++i)
        allRows.Add(i);
      double grandMean = this.TreeTargetMean(allRows);

      Node root = new Node();
      root.id = 0;
      root.value = grandMean;
      root.isLeaf = false;
      root.rows = allRows;
      this.tree[0] = root;

      // build tree breadth-first
      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node currNode = this.tree[i];
        if (currNode == null) continue;

        // safety checks
        if (currNode.id "gte" maxStartID ||
          currNode.rows.Count "lt" this.minSamples)
        {
          currNode.isLeaf = true;
          currNode.left = -1;   // Explicitly isolate node boundaries
          currNode.right = -1;
          currNode.colIdx = -1;
          continue;
        }

        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1];

        // check for split failure
        if (colIdx == -1)
        {
          currNode.isLeaf = true;
          currNode.left = -1;
          currNode.right = -1;
          currNode.colIdx = -1;
          continue;
        }

        // got a valid split point
        currNode.colIdx = colIdx;
        currNode.thresh = splitVal;

        // avoid continuous resizing allocations
        List"lt"int"gt" leftIdxs = 
          new List"lt"int"gt"(currNode.rows.Count);
        List"lt"int"gt" rightIdxs = 
          new List"lt"int"gt"(currNode.rows.Count);

        for (int k = 0; k "lt" currNode.rows.Count; ++k)
        {
          int r = currNode.rows[k];
          if (this.trainX[r][colIdx] "lte" splitVal)
            leftIdxs.Add(r);
          else
            rightIdxs.Add(r);
        }

        int leftID = currNode.id * 2 + 1;
        int rightID = currNode.id * 2 + 2;

        // check both branches
        bool leftValid = (leftID "lte" maxID && 
          leftIdxs.Count "gte" this.minLeaf);
        bool rightValid = (rightID "lte" maxID &&
          rightIdxs.Count "gte" this.minLeaf);

        if (leftValid == true && rightValid == true)
        {
          // create child nodes
          currNode.left = leftID;
          Node leftNode = new Node();
          leftNode.id = leftID;
          leftNode.rows = leftIdxs;
          leftNode.value = 
            this.TreeTargetMean(leftNode.rows);
          this.tree[leftID] = leftNode;

          currNode.right = rightID;
          Node rightNode = new Node();
          rightNode.id = rightID;
          rightNode.rows = rightIdxs;
          rightNode.value = 
            this.TreeTargetMean(rightNode.rows);
          this.tree[rightID] = rightNode;
        }
        else
        {
          // structural asymmetry/failure edge case
          // make parent into a leaf node
          currNode.isLeaf = true;
          currNode.left = -1;
          currNode.right = -1;
          currNode.colIdx = -1;
        }
      }

      // delete row info (for ensembles)
      if (this.saveRows == false)
      {
        for (int i = 0; i "lt" this.tree.Count; ++i)
          if (this.tree[i] != null)
            this.tree[i].rows = null;
      }
    }

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      int p = 0;
      double lastValidValue = 0.0;
      while (p != -1 && p "lt" this.tree.Count)
      {
        Node currNode = this.tree[p];

        if (currNode == null) break;

        lastValidValue = currNode.value;

        if (currNode.isLeaf == true ||
          currNode.colIdx == -1 ||
          currNode.left "gte" this.tree.Count ||
          currNode.right "gte" this.tree.Count)
          break;

        if (x[currNode.colIdx] "lte" currNode.thresh)
          p = currNode.left;
        else
          p = currNode.right;
      }
      return lastValidValue;
    }

    // ------------------------------------------------------

    public void Explain(double[] x)
    {
      int p = 0;
      double lastValidValue = 0.0;
      Node currNode = null;
      string s = "\nIF \n";

      while (p != -1 && p "lt" this.tree.Count)
      {
        currNode = this.tree[p];
        if (currNode == null) break;

        lastValidValue = currNode.value;
        if (currNode.isLeaf == true) break;

        s += "column " + currNode.colIdx + " ";

        if (x[currNode.colIdx] "lte" currNode.thresh)
        {
          s += " "lte" " +
            currNode.thresh.ToString("F4").PadLeft(8) +
            " AND \n";
          p = currNode.left;
        }
        else
        {
          s += " "gt"  " +
            currNode.thresh.ToString("F4").PadLeft(8) +
            " AND \n";
          p = currNode.right;
        }
      }

      int nid;
      if (currNode == null)
        nid = -1;
      else
        nid = currNode.id;

      s += "THEN node [" + nid + "] predicted = " +
        currNode.value.ToString("F4");
      Console.WriteLine(s);
    }

    // ------------------------------------------------------

    public double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" dataX.Length; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        if (Math.Abs(predY - actualY) "lt"
          (pctClose * Math.Abs(actualY)))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

    // ------------------------------------------------------

    public double MSE(double[][] dataX, double[] dataY)
    {
      // standard machine learning MSE, not tree MSE
      int n = dataX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        sum += (actualY - predY) * (actualY - predY);
      }
      return sum / n;
    }

    // ------------------------------------------------------

    public void Display()
    {
      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node n = this.tree[i];

        // check for empty nodes
        if (n == null) continue;

        string s1 = "ID " +
          n.id.ToString().PadRight(3) + " | ";
        string s2 = "sc " +
          n.colIdx.ToString().PadLeft(3) + " | ";
        string s3 = "sv " +
          n.thresh.ToString("F4").PadLeft(8) + " | ";
        string s4 = "L " +
          n.left.ToString().PadLeft(3) + " | ";
        string s5 = "R " +
          n.right.ToString().PadLeft(3) + " | ";
        string s6 = "py " +
          n.value.ToString("F4").PadLeft(8) + " | ";
        string s7 = "leaf " +
          (n.isLeaf == true ? "T" : "F") + " | ";
        string s8 = "rc " + n.rows.Count;

        Console.WriteLine(s1 + s2 + s3 + s4 +
          s5 + s6 + s7 + s8);
      }
    } // Display()

    // ------------------------------------------------------

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      // optimized for performance at expense of clarity
      int bestColIdx = -1;
      double bestThresh = 0.0;
      double bestVar = double.MaxValue;
      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

      if (nRows == 0)
      {
        throw new Exception("empty data in BestSplit()");
      }

      int[] colIndices = new int[nCols];
      for (int k = 0; k "lt" nCols; ++k)
        colIndices[k] = k;

      for (int i = 0; i "lt" nCols; ++i) // Fisher-Yates
      {
        int ri = rnd.Next(i, nCols);
        int tmp = colIndices[i];
        colIndices[i] = colIndices[ri];
        colIndices[ri] = tmp;
      }

      int numColsToUse;
      if (this.numSplitCols == -1)
        numColsToUse = nCols;
      else
        numColsToUse = Math.Min(this.numSplitCols, nCols);

      int[] activeCols = new int[numColsToUse];
      for (int k = 0; k "lt" numColsToUse; ++k)
        activeCols[k] = colIndices[k];

      int[] sortedRows = new int[nRows];
      double[] featureKeys = new double[nRows];

      double totalSum = 0.0;
      double totalSumSq = 0.0;
      for (int i = 0; i "lt" nRows; ++i)
      {
        double yCurr = this.trainY[rows[i]];
        totalSum += yCurr;
        totalSumSq += yCurr * yCurr;
      }

      for (int j = 0; j "lt" activeCols.Length; ++j)
      {
        int colIdx = activeCols[j];
        for (int i = 0; i "lt" nRows; ++i)
        {
          int r = rows[i];
          sortedRows[i] = r;
          featureKeys[i] = this.trainX[r][colIdx];
        }

        // sort both arrays based on keys
        Array.Sort(featureKeys, sortedRows);
        double leftSum = 0.0;
        double leftSumSq = 0.0;

        for (int i = 0; i "lt" nRows - 1; ++i)
        {
          int currRowIdx = sortedRows[i];
          double yCurr = this.trainY[currRowIdx];
          leftSum += yCurr;
          leftSumSq += yCurr * yCurr;

          int leftCount = i + 1;
          int rightCount = nRows - leftCount;

          double currFeatureVal = featureKeys[i];
          double nextFeatureVal = featureKeys[i + 1];

          // because x-features are sorted, this
          // code skips over previously seen values
          if (currFeatureVal == nextFeatureVal)
            continue;

          if (leftCount "lt" this.minLeaf ||
            rightCount "lt" this.minLeaf)
            continue;

          // tricky math to compute variances
          double rightSum = totalSum - leftSum;
          double rightSumSq = totalSumSq - leftSumSq;

          double tmp1 =
            (leftSum / leftCount) * (leftSum / leftCount);
          double leftVar = (leftSumSq / leftCount) - tmp1;

          double tmp2 =
            (rightSum / rightCount) * (rightSum / rightCount);
          double rightVar = (rightSumSq / rightCount) - tmp2;

          if (leftVar "lt" 0.0) leftVar = 0.0;
          if (rightVar "lt" 0.0) rightVar = 0.0;

          double weightedVar =
            ((leftCount * leftVar) + (rightCount * rightVar))
            / nRows;

          if (weightedVar "lt" bestVar)
          {
            bestColIdx = colIdx;
            // bestThresh =
            //   (currFeatureVal + nextFeatureVal) / 2.0;
            bestThresh = currFeatureVal; // simpler approach
            bestVar = weightedVar;
          }
        }
      }
      // using a tuple would create version dependency
      double[] result = new double[2];
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;
    }

    // ------------------------------------------------------

    private double TreeTargetMean(List"lt"int"gt" rows)
    {
      if (rows == null || rows.Count == 0) return 0.0;
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
      {
        sum += this.trainY[rows[i]];
      }
      return sum / rows.Count;
    }

  } // class DecisionTreeRegressor

} // ns

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

I ran into an interesting problem recently. I had a set of training data and I needed to condense it to a smaller size (number of items). This is a rare scenario because in most cases you want as much training data as possible.

There are many possible ways to prune a set of training data. There is no one best technique, and different techniques will give slightly different results. I decided to use an RBF kernel function. RBF (radial basis function) computes a measure of similarity between two vectors. RBF(x1, x2, gamma) gives a value between 0 (infinitely different vectors) and 1 (identical vectors). The gamma is a free parameter, typically around 0.5 or so. It must be tuned using trial and error.

I originally implemented a demo, with two different algorithms, using Python. The “KRA” (kernel matrix row average) is best understood by looking at the code. KRA is simple and deterministic, but relatively slow to compute. The “FFT” (farthest first traversal) is a standard algorithm. FFT is more complicated than KRA and is also non-deterministic, but is relatively faster to compute.

It was a bit trickier than I expected, but I eventually got the Python demo up and running. However, the project I’m working uses the C# language, so I wanted to refactor the Python pruning code to C#. I did so, but because Python/NumPy has hundreds of built-in functions and syntax such as np.argmin() and X_new = X[mask], I needed to implement a lot of C# utility code — about 150 lines.

After a couple of hours I had a C# version running. Sample output:

Begin prune training data using C# demo

Setting up dummy data

Source data:
   1.7641   0.4002   0.9787   2.2409
   1.8676  -0.9773   0.9501  -0.1514
  -0.1032   0.4106   0.1440   1.4543
   0.7610   0.1217   0.4439   0.3337
   1.4941  -0.2052   0.3131  -0.8541
  -2.5530   0.6536   0.8644  -0.7422
   2.2698  -1.4544   0.0458  -0.1872
   1.5328   1.4694   0.1549   0.3782
  -0.8878  -1.9808  -0.3479   0.1563
   9.0000   9.0000   9.0000   9.0000

Pruning down to 3 dissimilar items using KRA with
 RBF gamma = 0.5000
Done

Prune mask:
  9  5  8

Pruned dataset:
   9.0000   9.0000   9.0000   9.0000
  -2.5530   0.6536   0.8644  -0.7422
  -0.8878  -1.9808  -0.3479   0.1563

==========================

Pruning down to 3 dissimilar items using FFT with
 RBF gamma = 0.5000
Done

Prune mask:
  7  9  5

Pruned dataset:
   1.5328   1.4694   0.1549   0.3782
   9.0000   9.0000   9.0000   9.0000
  -2.5530   0.6536   0.8644  -0.7422

End demo

The demo data was generated randomly by the Python version. I manually modified the last row to make it wildly anomalous. The results of the C# pruning program were identical to the results of the Python program.

British secret agent James Bond pruned away many villains and villainesses. The first five Bond films, starring Sean Connery, had a profound effect on popular culture. “Dr. No” (1962), “From Russia with Love” (1963), “Goldfinger” (1964), “Thunderball” (1965), and “You Only Live Twice” (1967).

Right: In “Thunderball”, Fiona Volpe is an agent of SPECTRE who uses her beauty to ensnare and kill a NATO jet bomber pilot in order to steal two atomic bombs. She and Bond dance and she tries to position Bond to be shot by an assassin, but Bond sees the gunman, and he spins and Volpe is killed instead.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

namespace PruneTrainingDataset
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin prune training data" +
        " using C# demo ");

      Console.WriteLine("\nSetting up dummy data ");
      double[][] trainX = new double[10][];
      trainX[0] = 
        new double[] { 1.7641, 0.4002, 0.9787, 2.2409 };
      trainX[1] = 
        new double[] { 1.8676, -0.9773, 0.9501, -0.1514 };
      trainX[2] = 
        new double[] { -0.1032, 0.4106, 0.1440, 1.4543 };
      trainX[3] = 
        new double[] { 0.7610, 0.1217, 0.4439, 0.3337 };
      trainX[4] = 
        new double[] { 1.4941, -0.2052, 0.3131, -0.8541 };
      trainX[5] = 
        new double[] { -2.5530, 0.6536, 0.8644, -0.7422 };
      trainX[6] = 
        new double[] { 2.2698, -1.4544, 0.0458, -0.1872 };
      trainX[7] = 
        new double[] { 1.5328, 1.4694, 0.1549, 0.3782 };
      trainX[8] = 
        new double[] { -0.8878, -1.9808, -0.3479, 0.1563 };
      //trainX[9] = 
      //  new double[] { 1.2303, 1.2024, -0.3873, -0.3023 };
      trainX[9] = new double[] { 9, 9, 9, 9 }; // anomaly

      Console.WriteLine("\nSource data: ");
      Utils.MatShow(trainX, 4, 9);

      double gamma = 0.5;  // arbitrary RBF parameter
      
      Console.WriteLine("\nPruning down to 3 dissimilar " +
        "items using KRA with RBF gamma = " + 
        gamma.ToString("F4"));

      int[] pruneMask = Prune(trainX, 3, gamma);
      Console.WriteLine("Done ");

      Console.WriteLine("\nPrune mask: ");
      Utils.VecShow(pruneMask, 3);

      double[][] prunedX = 
        Utils.MatSelectRows(trainX, pruneMask);
      Console.WriteLine("\nPruned dataset: ");
      Utils.MatShow(prunedX, 4, 9);

      Console.WriteLine("\n========================== ");

      Console.WriteLine("\nPruning down to 3 dissimilar " +
        "items using FFT with RBF gamma = " + 
        gamma.ToString("F4"));

      Random rnd = new Random(0);
      pruneMask = PruneFFT(trainX, 3, gamma, rnd);
      Console.WriteLine("Done ");

      Console.WriteLine("\nPrune mask: ");
      Utils.VecShow(pruneMask, 3);

      prunedX = 
        Utils.MatSelectRows(trainX, pruneMask);
      Console.WriteLine("\nPruned dataset: ");
      Utils.MatShow(prunedX, 4, 9);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main

    // ------------------------------------------------------

    static int[] Prune(double[][] X, int nKeep,
      double gamma)
    {
      // Kernel matrix row average technique
      int n = X.Length;
      if (nKeep "gte" n)
        return Utils.VecRange(n);

      double[][] K = Utils.MatMake(n, n);  // all sims
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = i; j "lt" n; ++j)
        {
          double z = RBF(X[i], X[j], gamma);
          K[i][j] = z;
          K[j][i] = z;
        }
      }

      double[] rowSums = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        double currRowSum = 0.0;
        for (int j = 0; j "lt" n; ++j)
          currRowSum += K[i][j];
        rowSums[i] = currRowSum / n;
      }

      int[] sorted = Utils.ArgSort(rowSums);
      // extract first values
      int[] result = new int[nKeep];
      for (int i = 0; i "lt" nKeep; ++i)
        result[i] = sorted[i];
      return result;
    }

    // ------------------------------------------------------

    static int[] PruneFFT(double[][] X, int nKeep,
      double gamma, Random rnd)
    {
      // farthest first traversal (FFT) algorithm
      // faster than KRA but non-deterministic
      // use for moderate to large datasets

      int n = X.Length;
      if (nKeep "gte" n)
        return Utils.VecRange(n);

      // pick random item to start
      int firstIdx = rnd.Next(0, n);
      List"lt"int"gt" selectedIdxs = new List"lt"int"gt"();
      selectedIdxs.Add(firstIdx);

      // ----------------------------------------------------
      // local helper
      // ----------------------------------------------------

      double[] rbfSims(double[] vec, double[][] M,
        double gamma)
      {
        int nr = M.Length;
        int nc = M[0].Length;
        double[] result = new double[nr];
        for (int i = 0; i "lt" nr; ++i)
        {
          double sum = 0.0;
          for (int j = 0; j "lt" nc; ++j)
          {
            double diff = M[i][j] - vec[j];
            sum += diff * diff;
          }
          result[i] = Math.Exp(-1 * gamma * sum);
        }
        return result;
      }

      // ----------------------------------------------------

      double[] maxSims = rbfSims(X[firstIdx], X, gamma);
      while (selectedIdxs.Count "lt" nKeep)
      {
        int nextIdx = Utils.VecArgMin(maxSims);
        selectedIdxs.Add(nextIdx);
        double[] newSims = rbfSims(X[nextIdx], X, gamma);
        maxSims = Utils.VecMaximums(maxSims, newSims);
      }

      int[] result = selectedIdxs.ToArray();
      return result;
    }

    // ------------------------------------------------------

    // ------------------------------------------------------

    //static int[] Prune(double[][] X, int nKeep,
    //  double gamma, Random rnd)
    //{
    //  // farthest first traversal (FFT) algorithm


    //}

    // ------------------------------------------------------

    static double RBF(double[] v1, double[] v2,
      double gamma)
    {
      int n = v1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double d = v1[i] - v2[i];
        sum += d * d;
      }
      double result = Math.Exp(-1 * gamma * sum);
      return result;
    }

    // ------------------------------------------------------

  } // class Program

  // ========================================================

  public class Utils
  {
    // ------------------------------------------------------
    //
    // lots of helpers
    //
    // ------------------------------------------------------

    public static double[][] MatLoad(string fn,
      int[] usecols, char sep, string comment)
    {
      List"lt"double[]"gt" result = 
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    // ------------------------------------------------------

    public static double[] MatToVec(double[][] X)
    {
      int nRows = X.Length;
      int nCols = X[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = X[i][j];
      return result;
    }

    // ------------------------------------------------------

    public static void MatShow(double[][] m, int dec,
      int wid)
    {
      int nRows = m.Length; int nCols = m[0].Length;
      double small = 1.0 / Math.Pow(10, dec);
      for (int i = 0; i "lt" nRows; ++i)
      {
        for (int j = 0; j "lt" nCols; ++j)
        {
          double v = m[i][j];
          if (Math.Abs(v) "lt" small) v = 0.0;
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

    // ------------------------------------------------------

    public static double[][] MatSelectRows(double[][] X,
      int[] rows)
    {
      int nRowsSrc = X.Length;
      int nColsSrc = X[0].Length;
      int n = rows.Length;
      double[][] result = MatMake(n, nColsSrc);

      for (int i = 0; i "lt" n; ++i) // i pts into result
      {
        int srcRow = rows[i];
        for (int j = 0; j "lt" nColsSrc; ++j)
        {
          result[i][j] = X[srcRow][j];
        }
      }
      return result;
    }

    // ------------------------------------------------------

    public static double[][] MatSelectRows(double[][] X,
      List"lt"int"gt" rows)
    {
      return MatSelectRows(X, rows.ToArray());
    }

    // ------------------------------------------------------

    public static double[][] MatMake(int nRows, int nCols)
    {
      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols];
      return result;
    }

    // ------------------------------------------------------

    public static int[] VecRange(int n)
    {
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = i;
      return result;
    }

    // ------------------------------------------------------

    public static int[] ArgSort(double[] vec)
    {
      // doesn't modify vec
      int n = vec.Length;
      int[] idxs = new int[n];
      for (int i = 0; i "lt" n; ++i)
        idxs[i] = i;

      double[] dup = new double[n];
      for (int i = 0; i "lt" n; ++i)
        dup[i] = vec[i];

      Array.Sort(dup, idxs);  // sort idxs based on dup vals
      return idxs;
    }

    // ------------------------------------------------------

    public static int VecArgMin(double[] vec)
    {
      int minIdx = 0;
      double minVal = vec[0];
      for (int i = 0; i "lt" vec.Length; ++i)
      {
        if (vec[i] "lt" minVal)
        {
          minVal = vec[i];
          minIdx = i;
        }
      }
      return minIdx;
    }

    // ------------------------------------------------------

    public static double[] VecMaximums(double[] v1,
      double[] v2)
    {
      int n = v1.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        if (v1[i] "gt" v2[i])
          result[i] = v1[i];
        else
          result[i] = v2[i];
      }
      return result;
    }

    // ------------------------------------------------------

    public static void VecShow(double[] vec, 
      int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

    // ------------------------------------------------------

    public static void VecShow(int[] vec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString().PadLeft(wid));
      Console.WriteLine("");
    }

    // ------------------------------------------------------

  } // class Utils

} // ns

One Sunday morning, I figured I’d put together a demo four different different regularization techniques for linear regression: L1, L2, weight decay, and input noise. I have implemented all these techniques before, but never in the same program. I decided to use from-scratch Python.

For technical reasons, the four regularization techniques are best used in conjunction with stochastic gradient descent, rather than closed form MP pseudo-inverse training or closed form left pseudo-inverse via normal equations.

It’s important to note that L2 regularization, weight decay regularization, and input noise regularization are all mathematically equivalent, even though they are implemented differently.

Here are the key parts of the output of the demo:

Scratch linear regression with SGD training

Loading synthetic train (200) and test (40) data
Done

First three train X:
[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]
[ 0.0776 -0.1616  0.3704 -0.5911  0.7562]
[-0.9452  0.3409 -0.1654  0.1174 -0.7192]

First three train y:
0.4840
0.1568
0.8054

Using scikit LinearRegression module
Model weights:
[-0.2656  0.0333 -0.0454  0.0358 -0.1146]
Model bias = 0.3619
Accuracy train (within 0.10) = 0.4600
Accuracy test (within 0.10) = 0.6500
MSE train = 0.0026
MSE test = 0.0020

Basic SGD training (no regularization)
Early exit at epoch 408
Model weights:
[-0.2657  0.0335 -0.0455  0.0358 -0.1150]
Model bias = 0.3628
Accuracy train (within 0.10) = 0.4650
Accuracy test (within 0.10) = 0.6250
MSE train = 0.0026
MSE test = 0.0020

SGD with L2 training
Early exit at epoch 25
Model weights:
[-0.2646  0.0341 -0.0461  0.0369 -0.1144]
Model bias = 0.3614
Accuracy train (within 0.10) = 0.4600
Accuracy test (within 0.10) = 0.6750
MSE train = 0.0026
MSE test = 0.0020

SGD with L1 training
Model weights:
[-0.2633  0.0303 -0.0427  0.0327 -0.1124]
Model bias = 0.3646
Accuracy train (within 0.10) = 0.4900
Accuracy test (within 0.10) = 0.6500
MSE train = 0.0026
MSE test = 0.0020

SGD with weight decay training
Model weights:
[-0.2650  0.0327 -0.0448  0.0356 -0.1145]
Model bias = 0.3647
Accuracy train (within 0.10) = 0.4700
Accuracy test (within 0.10) = 0.6750
MSE train = 0.0026
MSE test = 0.0020

SGD with input noise training
Model weights:
[-0.2656  0.0334 -0.0470  0.0361 -0.1143]
Model bias = 0.3578
Accuracy train (within 0.10) = 0.4650
Accuracy test (within 0.10) = 0.7000
MSE train = 0.0026
MSE test = 0.0020

End demo

If you look closely, you’ll see that the MSE values for the four regularization techniques, and the two baseline techniques with no regularization, are all the same. This illustrates one of the main reasons why I almost never use regularization for linear regression with SGD training.

The key weight update statements for no regularization, L2 regularization, L1 regularization, weight decay regularization, and input noise regularization:

# 1. no regularization:
for j in range(dim):
  self.weights[j] -= lrn_rate * error * x[j]

# 2. L2:
for j in range(dim):
  self.weights[j] -= lrn_rate * error * x[j] + \
   (l2_lamda * self.weights[j])

# 3. L1:
for j in range(dim):
  self.weights[j] -= lrn_rate * error * x[j] + \
   (l1_lamda * np.sign(self.weights[j]))

# 4. weight decay:
for j in range(dim):
  self.weights[j] *= (1.0 - decay)
  self.weights[j] -= lrn_rate * error * x[j]

# 5. input noise:
x = X[idx]
x += self.rnd.normal(0.0, noise, dim)
. . . 
for j in range(dim):
  self.weights[j] -= lrn_rate * error * x[j]

From a theoretical point of view, regularization makes sense in order to keep the magnitude of the weight values from exploding. But from a practical point of view, regularization adds an extra hyperparameter to tune, and is often a waste of time.

If your goal is prediction on real-life data, using linear regression doesn’t make sense because almost no real-world data fits a purely linear model. And adding regularization isn’t going to improve the prediction model in any significant way.

If your goal is to use linear regression to establish a baseline result for comparison against more powerful regression techniques (quadratic regression, kernel ridge regression, neural network regression, gradient boost regression), then adding regularization is actually counter-productive — you want a basic linear model, not contaminated by regularization.

A fun mental exercise for me on a Sunday morning.

Linear regression is too weak for most practical scenarios, but it’s useful to establish baseline results to compare with more powerful techniques such as kernel (not colonel) ridge regression.

Right: In “The Invisible Boy” (1957), the plot is unbelievably strange, but briefly a supercomputer and a robot from the future go rogue. The computer is defeated and the robot turns good. Colonel Mackin is part of the unsuccessful military attempt to stop the supercomputer. This is a bad movie (my grade = C-), but the robot is the one used in “Forbidden Planet” (1956), one of the best science fiction movies of all time.

Demo program. Replace the instances of “lt” with the less-than Boolean operator symbol (my blod editor chokes on symbols).

# linear_regression_sgd_regularization_techniques.py

# basic SGD, L1 SGD, L2 SGD, wt-decay SGD, input-noise SGD
# note: L2, weight decay, and adding noise to input are
# all mathematically equivalent!

import numpy as np

# -----------------------------------------------------------

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed', linewidth=120)

# -----------------------------------------------------------
# external eval functions: accuracy(), mse().
# model has internal r2_score() method.
# -----------------------------------------------------------

def accuracy(model, data_X, data_y, pct_close):
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    if np.abs(y - pred_y) "lt" np.abs(y * pct_close):
      n_correct += 1
    else: 
      n_wrong += 1
  return n_correct / (n_correct + n_wrong)

# -----------------------------------------------------------

def mse(model, data_X, data_y):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    diff = pred_y - y
    sum += diff * diff
  return sum /n

# ===========================================================

class LinearRegressionSGD:
  def __init__(self, seed=0):
    self.weights = None
    self.bias = None
    self.rnd = np.random.RandomState(seed)

  # ---------------------------------------------------------

  def fit_basic(self, X, y, lrn_rate, max_epochs, exit_tol):
    # no regulariation
    n = len(X); dim = len(X[0])
    self.weights = np.zeros(dim)
    self.bias = 0.0

    indices = np.arange(n)
    consecutive_no_change = 0
    for epoch in range(max_epochs):
      self.rnd.shuffle(indices)
      weights_old = self.weights.copy()

      for idx in indices:
        x = X[idx]
        target_y = y[idx]

        pred_y = 0.0
        for j in range(dim):
          pred_y += self.weights[j] * x[j]
        pred_y += self.bias

        error = pred_y - target_y
        for j in range(dim):
          self.weights[j] -= lrn_rate * error * x[j]
        self.bias -= lrn_rate * error

      wts_change = \
        self.euc_distance(self.weights, weights_old)
      if wts_change "lt" exit_tol:
        consecutive_no_change += 1
        if consecutive_no_change == 3:
          print("Early exit at epoch " + str(epoch))
          break
      else:
        consecutive_no_change = 0

  # ---------------------------------------------------------

  def fit_L2(self, X, y, lrn_rate, max_epochs, \
    l2_lamda, exit_tol):
    n = len(X); dim = len(X[0])
    self.weights = np.zeros(dim)
    self.bias = 0.0

    indices = np.arange(n)
    for epoch in range(max_epochs):
      self.rnd.shuffle(indices)
      weights_old = self.weights.copy()

      for idx in indices:
        x = X[idx]
        target_y = y[idx]

        pred_y = 0.0
        for j in range(dim):
          pred_y += self.weights[j] * x[j]
        pred_y += self.bias

        error = pred_y - target_y
        for j in range(dim):
          self.weights[j] -= lrn_rate * error * x[j] + \
            (l2_lamda * self.weights[j])
        self.bias -= lrn_rate * error

      wts_change = \
        self.euc_distance(self.weights, weights_old)
      if wts_change "lt" exit_tol:
        print("Early exit at epoch " + str(epoch))
        break

  # ---------------------------------------------------------

  def fit_L1(self, X, y, lrn_rate, max_epochs, \
    l1_lamda, exit_tol):
    n = len(X); dim = len(X[0])
    self.weights = np.zeros(dim)
    self.bias = 0.0

    indices = np.arange(n)
    for epoch in range(max_epochs):
      self.rnd.shuffle(indices)
      weights_old = self.weights.copy()

      for idx in indices:
        x = X[idx]
        target_y = y[idx]

        pred_y = 0.0
        for j in range(dim):
          pred_y += self.weights[j] * x[j]
        pred_y += self.bias

        error = pred_y - target_y
        for j in range(dim):
          self.weights[j] -= lrn_rate * error * x[j] + \
            (l1_lamda * np.sign(self.weights[j]))
        self.bias -= lrn_rate * error

      wts_change = \
        self.euc_distance(self.weights, weights_old)
      if wts_change "lt" exit_tol:
        print("Early exit at epoch " + str(epoch))
        break

  # ---------------------------------------------------------

  def fit_decay(self, X, y, lrn_rate, max_epochs, \
    decay, exit_tol):
    n = len(X); dim = len(X[0])
    self.weights = np.zeros(dim)
    self.bias = 0.0

    indices = np.arange(n)
    for epoch in range(max_epochs):
      self.rnd.shuffle(indices)
      weights_old = self.weights.copy()

      for idx in indices:
        x = X[idx]
        target_y = y[idx]

        pred_y = 0.0
        for j in range(dim):
          pred_y += self.weights[j] * x[j]
        pred_y += self.bias

        error = pred_y - target_y
        for j in range(dim):
          self.weights[j] *= (1.0 - decay)  # note
          self.weights[j] -= lrn_rate * error * x[j]
        self.bias -= lrn_rate * error

      wts_change = \
        self.euc_distance(self.weights, weights_old)
      if wts_change "lt" exit_tol:
        print("Early exit at epoch " + str(epoch))
        break

  # ---------------------------------------------------------

  def fit_input_noise(self, X, y, lrn_rate, max_epochs, \
    noise, exit_tol):
    n = len(X); dim = len(X[0])
    self.weights = np.zeros(dim)
    self.bias = 0.0

    indices = np.arange(n)
    for epoch in range(max_epochs):
      self.rnd.shuffle(indices)
      weights_old = self.weights.copy()

      for idx in indices:
        x = X[idx]
        x += self.rnd.normal(0.0, noise, dim) # note
        target_y = y[idx]

        pred_y = 0.0
        for j in range(dim):
          pred_y += self.weights[j] * x[j]
        pred_y += self.bias

        error = pred_y - target_y
        for j in range(dim):
          self.weights[j] -= lrn_rate * error * x[j]
        self.bias -= lrn_rate * error

      wts_change = \
        self.euc_distance(self.weights, weights_old)
      if wts_change "lt" exit_tol:
        print("Early exit at epoch " + str(epoch))
        break

  # ---------------------------------------------------------

  def euc_distance(self, v1, v2):
    n = len(v1)
    sum = 0.0
    for i in range(n):
      sum += (v1[i] - v2[i]) * (v1[i] - v2[i])
    result = np.sqrt(sum)
    return result

  # ---------------------------------------------------------

  def predict_one(self, x):
    n = len(x)
    sum = 0.0
    for i in range(n):
      sum += self.weights[i] * x[i]
    sum += self.bias
    return sum

  # ---------------------------------------------------------

  def predict(self, X):
    n = len(X)
    result = np.zeros(n)
    for i in range(n):
      result[i] = self.predict_one(X[i])
    return result

# ===========================================================

def main():
  print("\nScratch linear regression with SGD training ")

  print("\nLoading synthetic train (200) and test (40) data")
  train_Xy = np.loadtxt(".\\Data\\synthetic_train_200.txt",
    usecols=[0,1,2,3,4,5], delimiter=",")
  train_X = train_Xy[:,[0,1,2,3,4]]
  train_y = train_Xy[:,5]

  test_Xy = np.loadtxt(".\\Data\\synthetic_test_40.txt",
    usecols=[0,1,2,3,4,5], delimiter=",")
  test_X = test_Xy[:,[0,1,2,3,4]]
  test_y = test_Xy[:,5]
  print("Done ")

  print("\nFirst three train X: ")
  for i in range(3):
    print(train_X[i])
  print("\nFirst three train y: ")
  for i in range(3):
    print("%0.4f " % train_y[i])

  # ---------------------------------------------------------

  from sklearn.linear_model import LinearRegression
  print("\nUsing scikit LinearRegression module ")
  model = LinearRegression()
  model.fit(train_X, train_y)
  print("Model weights: ")
  print(model.coef_)
  print("Model bias = %0.4f " % model.intercept_)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  mse_train = mse(model, train_X, train_y)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test)

  # ---------------------------------------------------------

  print("\nBasic SGD training (no regularization) ")
  model = LinearRegressionSGD(seed=0)
  lrn_rate = 0.01
  max_epochs = 1000
  exit_tol = 0.001

  model.fit_basic(train_X, train_y, lrn_rate, \
    max_epochs, exit_tol)

  print("Model weights: ")
  print(model.weights)
  print("Model bias = %0.4f " % model.bias)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  mse_train = mse(model, train_X, train_y)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test)

  # ---------------------------------------------------------

  print("\nSGD with L2 training ")
  model = LinearRegressionSGD(seed=0)
  lrn_rate = 0.01
  max_epochs = 1000
  lamda = 0.00001
  exit_tol = 0.001

  model.fit_L2(train_X, train_y, lrn_rate, \
    max_epochs, lamda, exit_tol)

  print("Model weights: ")
  print(model.weights)
  print("Model bias = %0.4f " % model.bias)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  mse_train = mse(model, train_X, train_y)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test) 

  # ---------------------------------------------------------

  print("\nSGD with L1 training ")
  model = LinearRegressionSGD(seed=0)
  lrn_rate = 0.01
  max_epochs = 1000
  lamda = 0.00001
  exit_tol = 0.0001

  model.fit_L1(train_X, train_y, lrn_rate, \
    max_epochs, lamda, exit_tol)

  print("Model weights: ")
  print(model.weights)
  print("Model bias = %0.4f " % model.bias)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  mse_train = mse(model, train_X, train_y)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test) 

  # ---------------------------------------------------------

  print("\nSGD with weight decay training ")
  model = LinearRegressionSGD(seed=0)
  lrn_rate = 0.01
  max_epochs = 1000
  decay = 0.00001
  exit_tol = 0.0001

  model.fit_decay(train_X, train_y, lrn_rate, \
    max_epochs, decay, exit_tol)

  print("Model weights: ")
  print(model.weights)
  print("Model bias = %0.4f " % model.bias)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  mse_train = mse(model, train_X, train_y)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test) 

  # ---------------------------------------------------------

  print("\nSGD with input noise training ")
  model = LinearRegressionSGD(seed=0)
  lrn_rate = 0.01
  max_epochs = 1000
  noise = 0.0001
  exit_tol = 0.0001

  model.fit_input_noise(train_X, train_y, lrn_rate, \
    max_epochs, noise, exit_tol)

  print("Model weights: ")
  print(model.weights)
  print("Model bias = %0.4f " % model.bias)
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("Accuracy train (within 0.10) = %0.4f " % acc_train)
  mse_train = mse(model, train_X, train_y)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Accuracy test (within 0.10) = %0.4f " % acc_test)
  print("MSE train = %0.4f " % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("MSE test = %0.4f " % mse_test) 

  # ---------------------------------------------------------

  print("\nEnd demo ")

# -----------------------------------------------------------

if __name__ == "__main__":
  main()

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

One morning before work, I figured I’d implement kernel support vector regression (SVR), from scratch, using Python. It took me a bit longer than expected but I got a demo up and running.

My implementation uses hard-wired RBF (radial basis function) as the kernel function. I use the sequential minimal optimization (SMO) algorithm for training.

I won’t try to explain SVR or SMO, but I will caution you that I’ve seen a huge amount of grossly incorrect information about these two topics.

Here’s the output of my demo:

Begin scratch Python SVR using SMO training

Loading synthetic train (200) and test (40) data
Done

First three train X:
[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]
[ 0.0776 -0.1616  0.3704 -0.5911  0.7562]
[-0.9452  0.3409 -0.1654  0.1174 -0.7192]

First three train y:
0.4840
0.1568
0.8054

Creating SVR-SMO model
Setting gamma = 0.3000
Setting C = 1.00
Setting epsilon = 0.0325
Setting max_iter = 100
Setting KKT tol = 0.001000

Creating and training SVR model using SMO
Done

Model dual coefs:
[ 0.1209  0.2127  0.2883 -0.2028  . . .  0.0180
  0.1330  0.0164  0.0507 -0.1139  . . .  0.1208
  . . .
  0.0778  0.0860  0.0032  0.0283  . . .  0.0015
 -0.0874 -0.2291  0.0082 -0.0202  . . .  0.0391
  0.2399]

Model bias = 0.4123

Number support vectors = 197

Train accuracy (0.10) = 0.9200
Test accuracy (0.10) = 0.9250

Train MSE = 0.0002
Test MSE = 0.0002

Train R2 = 0.9951
Test R2 = 0.9941

Predicting for train_X[0]
Predicted y = 0.4913

End demo

The demo data is synthetic. It was generated by a 5-10-1 neural network with random weights and biases. There are 200 training items and 40 test items. The data is quasi-normalized, which is needed in most situations when using SVR.

I don’t like to be negative, but I am not a fan of SVR. Kernel ridge regression (KRR), which is closely related to SVR, is clearly superior to SVR. SVR is more complicated to implement than KRR. SVR is much more difficult to train than KRR. SVR usually gives less accuracy than KRR (mostly because of the difficulty tuning SVR parameters). But there are a few relatively rare scenarios where SVR works well.

I used the standard SVR SMO training algorithm design with two parallel alpha weight vectors, named alpha and alpha*. These are Lagrange multipliers associated with data points that fall above the upper margin of the epsilon error tube, and points that fall below the epsilon tube. It is possible to use a single alpha vector, but my experiments showed that such a design seems to always retain all training items as the support vectors.

The SMO algorithm is very complex, but the key parts look like:

init alpha[] and alpha*[] to number of training items
loop several times
  loop each train item i
    compute predicted y and error for i
    check "KKT" conditions to see if i is "good"
    if i is good
      pick a random item j that is not i
      compute predicted y and error for j
      use very complex math to update alpha[i], alpha[j],
        alpha*[i], alpha*[j]
    else
      next i
    end-if
  end-loop
end-loop

Because the SMO algorithm is so complex, all of the implementations I’ve seen are wildly different from each other.

Good fun.

Support vector regression made a big splash in the late 1990s, but faded out quickly as soon as people discovered all of its problems, and realized the clear superiority of the closely-related kernel ridge regression (KRR) for most problem scenarios. But I admire the passion that some people have for SVR. When I worked at Microsoft Research, my office in Building 99 was just a few doors away from the office of John Platt, the inventor of the SMO algorithm. Platt is at Google now.

Right: “Josh Kirby: Time Warrior” (1995) tells the story of . . . well, I really couldn’t quite follow it. Josh Kirby is a 14-year-old boy who is visited by people from the 25th century. Over the course of six 90-minute videos, they travel through time to locate the pieces of a super weapon. The stop motion effects are quite good . . . the movie, not so much. But again, I give the movie creators a lot of credit for their passion.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. My blog editor chokes on symbols.

# svr_smo.py

# kernel support vector regression from scratch.
# uses a hard-wired RBF kernel function.

import numpy as np

# -----------------------------------------------------------

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed', linewidth=120)

# -----------------------------------------------------------
# external eval functions: accuracy(), mse(), r2_score()
# -----------------------------------------------------------

def accuracy(model, data_X, data_y, pct_close):
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    if np.abs(y - pred_y) "lt" np.abs(y * pct_close):
      n_correct += 1
    else: 
      n_wrong += 1
  return n_correct / (n_correct + n_wrong)

# -----------------------------------------------------------

def mse(model, data_X, data_y):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    diff = pred_y - y
    sum += diff * diff
  return sum /n

# -----------------------------------------------------------

def r2_score(model, data_X, data_y):
  # coefficient of determination == scikit score()
  ss_res = 0.0
  ss_tot = 0.0
  n = len(data_X)
  mean_y = np.mean(data_y)
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    ss_res += (y - pred_y) * (y - pred_y)
    ss_tot += (y - mean_y) * (y - mean_y)
  result = 1.0 - (ss_res / ss_tot)
  return result

# ===========================================================

class KernelSVR:
  def __init__(self, gamma=0.1, epsilon=0.1, C=1.0,
    max_iter=100, tol=1.0e-3, seed=1):
    self.gamma = gamma
    self.epsilon = epsilon
    self.C = C
    self.max_iter = max_iter
    self.tol = tol  # for KKT conditions
    self.rnd = np.random.RandomState(seed)
    
    self.alpha = None
    self.alpha_star = None
    self.b = 0.0
    self.supp_X = None
    self.supp_y = None
    self.dual_weights = None
  # ---------------------------------------------------------

  def kernel_matrix(self, X1, X2):
    sq_dist = np.sum(X1**2, axis=1).reshape(-1, 1) + \
      np.sum(X2**2, axis=1) - 2 * np.dot(X1, X2.T)
    return np.exp(-self.gamma * sq_dist)

  # ---------------------------------------------------------

  def fit(self, X, y):
    n = X.shape[0]
    self.alpha = np.zeros(n)
    self.alpha_star = np.zeros(n)
    self.b = np.mean(y)
    
    K = self.kernel_matrix(X, X)
    n_passes = 0
    while n_passes "lt" self.max_iter:
      num_changed_alphas = 0
      
      for i in range(n):
        # prediction and error for i
        pred_i = np.dot(self.alpha_star - \
          self.alpha, K[:, i]) + self.b
        err_i = pred_i - y[i]
        
        # check KKT conditions for SVR within tolerance
        if ((err_i "gt" self.epsilon + self.tol and \
             self.alpha[i] "lt" self.C) or \
            (err_i "gt" self.epsilon + self.tol and \
             self.alpha_star[i] "gt" 0) or
            (err_i "lt" -self.epsilon - self.tol and \
             self.alpha[i] "gt" 0) or \
            (err_i "lt" -self.epsilon - self.tol and \
             self.alpha_star[i] "lt" self.C)):
            
          # pick a random second index j != i
          j = i
          while j == i:
            j = self.rnd.randint(0, n)
              
          pred_j = np.dot(self.alpha_star - \
            self.alpha, K[:, j]) + self.b
          err_j = pred_j - y[j]
          
          # save old values
          a_i_old, as_i_old = self.alpha[i], \
            self.alpha_star[i]
          a_j_old, as_j_old = self.alpha[j], \
            self.alpha_star[j]
          
          # kernel second derivative step denominator
          eta = K[i, i] + K[j, j] - 2 * K[i, j]
          if eta "lte" 0:
            continue
              
          # SVR linear constraint constant
          constraint = (as_i_old - a_i_old) + (as_j_old - a_j_old)
          
          # joint proxy variable s_j = (alpha_j* - alpha_j)
          s_j_old = as_j_old - a_j_old
          s_j_new = s_j_old + (err_i - err_j) / eta
          
          # bounds L and H for the proxy s_j
          L = max(-self.C, constraint - self.C)
          H = min(self.C, constraint + self.C)
          s_j_new = np.clip(s_j_new, L, H)
          
          if abs(s_j_new - s_j_old) "lt" 1.0e-5:
            continue
              
          # reconstruct alpha_j and alpha_j* 
          if s_j_new "gte" 0.0:
            self.alpha_star[j] = s_j_new
            self.alpha[j] = 0.0
          else:
            self.alpha_star[j] = 0.0
            self.alpha[j] = -s_j_new
              
          # update alpha_i and alpha_i*
          s_i_new = constraint - s_j_new
          if s_i_new "gte" 0.0:
            self.alpha_star[i] = s_i_new
            self.alpha[i] = 0.0
          else:
            self.alpha_star[i] = 0.0
            self.alpha[i] = -s_i_new
              
          # update bias
          b1 = self.b - err_i - ((self.alpha_star[i] - \
               self.alpha[i]) - (as_i_old - a_i_old)) * \
               K[i, i] - ((self.alpha_star[j] - \
               self.alpha[j]) - (as_j_old - a_j_old)) * \
               K[i, j]
          b2 = self.b - err_j - ((self.alpha_star[i] - \
               self.alpha[i]) - (as_i_old - a_i_old)) * \
               K[i, j] - ((self.alpha_star[j] - \
               self.alpha[j]) - (as_j_old - a_j_old)) * \
               K[j, j]
          
          if 0.0 "lt" self.alpha[i] "lt" self.C or \
            0.0 "lt" self.alpha_star[i] "lt" self.C:
            self.b = b1
          elif 0 "lt" self.alpha[j] "lt" self.C or \
            0.0 "lt" self.alpha_star[j] "lt" self.C:
            self.b = b2
          else:
            self.b = (b1 + b2) / 2.0
              
          num_changed_alphas += 1
          
      if num_changed_alphas == 0:
        n_passes += 1
      else:
        n_passes = 0
        
    # prune: store only explicit support vectors

    # 1. combine alpha and alpha*
    self.dual_weights = self.alpha_star - self.alpha

    # 2. compute mask
    sv_mask = np.abs(self.dual_weights) "gt" 1.0e-4

    # 3. mask supp vecs
    self.supp_X = X[sv_mask]
    self.supp_y = y[sv_mask]

    # 4. mask weights
    self.dual_weights = self.dual_weights[sv_mask]

    return self
  
  # ---------------------------------------------------------  

  def predict(self, X):
    K = self.kernel_matrix(X, self.supp_X)
    return np.dot(K, self.dual_weights) + self.b

  # ---------------------------------------------------------

  # def get_supp_idxs(self):
  #   TODO  result = []
  #   return result

# ===========================================================

def main():
  print("\nBegin scratch Python SVR using SMO training ")

  ## quick sanity check
  # np.random.seed(0)
  # n_samples = 40; n_features = 4
  # X = np.random.randn(n_samples, n_features)
  # weights = np.array([0.2, -0.5, 0.3,  0.1])
  # bias = 0.45
  # y = X @ weights + bias + np.random.randn(n_samples)

  # print("\nX = "); print(X[0:3,:]); print(" . . . ")
  # print("\ny = "); print(y[0:3], end=""); print(" . . . ")

  # model = KernelSVR(gamma=0.50, epsilon=0.01, C=1.0, 
  #   max_iter=20, tol=1.0e-5)
  # model.fit(X, y)

  # MSE = mse(model, X, y)
  # print("\nModel MSE = %0.4f " % MSE)

  print("\nLoading synthetic train (200) and test (40) data")
  train_Xy = np.loadtxt(".\\Data\\synthetic_train_200.txt",
    usecols=[0,1,2,3,4,5], delimiter=",")
  train_X = train_Xy[:,[0,1,2,3,4]]
  train_y = train_Xy[:,5]

  test_Xy = np.loadtxt(".\\Data\\synthetic_test_40.txt",
    usecols=[0,1,2,3,4,5], delimiter=",")
  test_X = test_Xy[:,[0,1,2,3,4]]
  test_y = test_Xy[:,5]
  print("Done ")

  print("\nFirst three train X: ")
  for i in range(3):
    print(train_X[i])
  print("\nFirst three train y: ")
  for i in range(3):
    print("%0.4f " % train_y[i])

  # ** SCIKIT results **
  # Setting gamma = 0.3000
  # Setting C = 1.0
  # Setting epsilon = 0.0010
  # Number model support vectors: [184]
  # Model bias: 0.4063
  # Train accuracy (0.10) = 0.9850
  # Test accuracy (0.10) = 0.9500  
  # Train MSE = 0.0000
  # Test MSE = 0.0002
  # Train R2 = 0.9988
  # Test R2 = 0.9930

  # create and train model
  print("\nCreating SVR-SMO model ")
  gamma = 0.20
  epsilon = 0.0325
  C = 1.0
  max_iter = 100  # max number iter with no improve
  tol = 1.0e-3

  print("Setting gamma = %0.4f " % gamma)
  print("Setting C = %0.2f " % C)
  print("Setting epsilon = %0.4f " % epsilon)
  print("Setting max_iter = " + str(max_iter))
  print("Setting tol = %0.6f " % tol)

  print("\nCreating and training SVR model using SMO ")

  model = KernelSVR(gamma=gamma, epsilon=epsilon, C=C, 
    max_iter=max_iter, tol=tol)

  model.fit(train_X, train_y)
  print("Done ")

  # print("\nModel alpha: ")
  # print(model.alpha)
  # print("\nModel alpha*: ") 
  # print(model.alpha_star)

  print("\nModel dual coefs: ")
  print(model.dual_weights)

  print("\nModel bias = %0.4f " % model.b)

  supp_vec_idxs = model.get_supp_idxs()
  print("\nNumber support vectors = " + \
    str(len(model.dual_weights)))

  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("\nTrain accuracy (0.10) = %0.4f" % acc_train)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Test accuracy (0.10) = %0.4f" % acc_test)

  mse_train = mse(model, train_X, train_y)
  print("\nTrain MSE = %0.4f" % mse_train)
  mse_test = mse(model, test_X, test_y)
  print("Test MSE = %0.4f" % mse_test)

  r2_train = r2_score(model, train_X, train_y)
  print("\nTrain R2 = %0.4f" % r2_train)
  r2_test = r2_score(model, test_X, test_y)
  print("Test R2 = %0.4f" % r2_test)

  print("\nPredicting for train_X[0] ")
  x = train_X[0].reshape(1,-1)
  pred_y = model.predict(x)[0]
  print("Predicted y = %0.4f " % pred_y)

  print("\nEnd demo ")

# -----------------------------------------------------------

if __name__ == "__main__":
  main()

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

I ran into an interesting problem recently. I had a set of training data and I needed to condense it to a smaller size (number of items). This is a rare scenario: in most cases you want as much training data as possible. For my scenario, I wanted to remove items that are very similar, or equivalently, retain items that are most dissimilar. Additionally, because the project I was working on used RBF similarity, I wanted to use RBF similarity to prune the dataset.

RBF (radial basis function) computes a measure of similarity between two vectors. RBF(x1, x2, gamma) gives a value between 0 (infinitely different) and 1 (the same). The gamma is a free parameter, typically around 0.5 or so.

There are many possible ways to prune a set of training data, and different algorithms will give slightly different results. There is no single best approach.

I came up with two different RBF-based pruning functions. The first one is principled, but complicated, and is non-deterministic. It is called the farthest first traversal (FFT) algorithm. In high-level pseudo-code:

prune_fft:
pick a random index, add it to list of selected items
for i = 1 to number items desired
  use RBF kernel function to find most dissimilar item to curr set
  add the dissimilar item to selected indexes
end-for
return list of selected indexes

The idea is a bit subtle. By adding the most dissimilar item to the current result set, you avoid adding similar items which are somewhat redundant.

The second pruning function is simpler, and is deterministic, but is slower than FFT. I call it kernel row average (KRA). In pseudo-code:

prune_kr1:
compute all pairs of RBF similarities
compute each row average, as similarity to all other items
sort row averages from low to high
extract first number items desired

So, row [0] of the RBF Kernel matrix holds similarity of item [0] to item [0], [1], [2] . . And row [1] of the RBF Kernel matrix holds similarity of item [1] to item [0], [1], [2] . . And so on. If you compute the average of each row, you get an average similarity for each item. If you sort those averages from low to high, the first n_to_keep items are the most dissimilar and the ones to retain.

I implemented a demo using Python. It was a bit trickier than I expected, but I eventually got the demo up and running. To torture myself, I made two versions of the FFT prune function, and two versions of the KRA prune function (one tricky but efficient, one clear but less efficient. Sample output.

Begin prune training data demo

Generating dummy training data
Done

Source data:
[[ 1.7641  0.4002  0.9787  2.2409]
 [ 1.8676 -0.9773  0.9501 -0.1514]
 [-0.1032  0.4106  0.1440  1.4543]
 [ 0.7610  0.1217  0.4439  0.3337]
 [ 1.4941 -0.2052  0.3131 -0.8541]
 [-2.5530  0.6536  0.8644 -0.7422]
 [ 2.2698 -1.4544  0.0458 -0.1872]
 [ 1.5328  1.4694  0.1549  0.3782]
 [-0.8878 -1.9808 -0.3479  0.1563]
 [ 9.0000  9.0000  9.0000  9.0000]]

Pruning down to 3 dissimilar items with RBF gamma = 0.5000

============================
1. Using FFT - tricky RBF
Prune mask:
[5, 9, 6]
============================
2. Using FFT - clear RBF
Prune mask:
[5, 9, 6]
============================
3. Using KRA - full memory K
Prune mask:
[9 5 8]
============================
4. Using KRA - low memory K
Prune mask:
[9 5 8]
============================

End demo

I manually inserted the last row with all 9.0 values to create one item that is clearly wildly different. The FFT and KRE pruning functions gave slightly different results, as expected. It was an interesting little exploration.

For my dataset pruning function, I used RBF kernel similarity, but I could have used Euclidean distance, or many other measures of similarity/dissimilarity.

Right: Mug shots of two random criminals. I could never tell them apart in a suspect lineup.

Demo program. Replace “lt” and “gte” with Boolean operator symbols.

# pruning_demo.py
# prune dataset to one with dissimilar items

import numpy as np

# -----------------------------------------------------------

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed', linewidth=120)

# -----------------------------------------------------------

def prune_FFT_1(X, n_to_keep, rnd, gamma=1.0):
  # returns indices into X of reduced version of X
  # where the items are most dissimilar so that nearly
  # duplicate rows are effectively removed.
  # use farthest first traversal (FFT) algorithm
  # fast but non-deterministic.
  # uses 'tricky' nested helper for RBF

  n = X.shape[0]
  if n_to_keep "gte" n:
    return np.arange(n)    
    
  # pick random item to start
  first_idx = rnd.randint(0, n)
  selected_idxs = [first_idx]  # a list

  # ---------------------------------------------------------
  # nested helper function
  # ---------------------------------------------------------

  def rbf_sims(row, M, gamma):
    # similarities between a row in M and all rows in M
    # efficient but tricky Python syntax
    # tmp: calculate squared Euclidean distances
    sq_dist = np.sum((M - row) ** 2, axis=1)
    result = np.exp(-1 * gamma * sq_dist)
    return result

  # ---------------------------------------------------------

  max_sims = rbf_sims(X[first_idx], X, gamma)
  while len(selected_idxs) "lt" n_to_keep:
    next_idx = np.argmin(max_sims).item()
    selected_idxs.append(next_idx)
    new_sims = rbf_sims(X[next_idx], X, gamma)
    max_sims = np.maximum(max_sims, new_sims)

  return selected_idxs  

# -----------------------------------------------------------

def prune_FFT_2(X, n_to_keep, rnd, gamma=1.0):
  # returns indices into X of reduced version of X
  # where the items are most dissimilar so that nearly
  # duplicate rows are effectively removed.
  # uses 'clear' nested helper

  n = X.shape[0]
  if n_to_keep "gte" n:
    return np.arange(n)    
    
  # pick random item to start
  first_idx = rnd.randint(0, n)
  selected_idxs = [first_idx]  # a list

  # ---------------------------------------------------------
  # nested helper function
  # ---------------------------------------------------------

  def rbf_sims(row, M, gamma):
    # similarities between a row in M and all rows in M
    # inefficient, but clear syntax
    n = len(M); dim = len(M[0])
    result = np.zeros(n)
    for i in range(n):
      sum = 0.0
      for j in range(dim):
        sum += (M[i][j] - row[j]) * (M[i][j] - row[j])
      result[i] = np.exp(-1 * gamma * sum)
    return result

  # -------------------------------------------------------

  max_sims = rbf_sims(X[first_idx], X, gamma)
  while len(selected_idxs) "lt" n_to_keep:
    next_idx = np.argmin(max_sims).item()
    selected_idxs.append(next_idx)
    new_sims = rbf_sims(X[next_idx], X, gamma)
    max_sims = np.maximum(max_sims, new_sims)

  return selected_idxs  

# -----------------------------------------------------------

def prune_KRA_1(X, n_to_keep, gamma=1.0):
  # use Kernel row averages. stores full K.
  # deterministic but slower than FFT
  n = len(X)
  if n_to_keep "gte" n:
    return np.arange(n)

  # ---------------------------------------------------------
  # nested helper
  # ---------------------------------------------------------

  def rbf(v1, v2, gamma):
    n = len(v1)
    sum =0.0
    for i in range(n):
      sum += (v1[i] - v2[i]) * (v1[i] - v2[i])
    return np.exp(-1 * gamma * sum)

  # ---------------------------------------------------------
  K = np.zeros((n,n))  # all similarity pairs
  for i in range(n):
    for j in range(i,n):
      z = rbf(X[i], X[j], gamma)
      K[i,j] = z; K[j,i] = z
  row_sums = np.zeros(n)
  for i in range(n):
    row_sum = 0.0
    for j in range(n):
      row_sum += K[i,j]
    row_sums[i] = row_sum / n

  sorted_sims = np.argsort(row_sums)  # small to large
  result = sorted_sims[0:n_to_keep] # first few
  return result

# -----------------------------------------------------------

def prune_KRA_2(X, n_to_keep, gamma=1.0):
  # this version doesn't store a large K matrix
  # but recomputes over and over
  n = len(X)
  if n_to_keep "gte" n:
    return np.arange(n)

  # ---------------------------------------------------------
  # nested helper
  # ---------------------------------------------------------

  def rbf(v1, v2, gamma):
    n = len(v1)
    sum =0.0
    for i in range(n):
      sum += (v1[i] - v2[i]) * (v1[i] - v2[i])
    return np.exp(-1 * gamma * sum)

  # ---------------------------------------------------------
  
  sims = np.zeros(n)
  for i in range(n):
    for j in range(n):
      sims[i] += rbf(X[i], X[j], gamma)
    sims[i] /= n
  sorted_sims = np.argsort(sims)  # small to large
  # small values are dissimilar
  result = sorted_sims[0:n_to_keep] # first values
  return result   

# -----------------------------------------------------------

print("\nBegin prune training data demo ")

rnd = np.random.RandomState(0)
print("\nGenerating dummy training data")
X = rnd.randn(10, 4)  # 10-by-3
print("Done ")
X[9] = np.array([9,9,9,9]) # make wildly different

print("\nSource data: ")
print(X)

gamma = 0.5
n_to_keep = 3

print("\nPruning down to 3 dissimilar items " + \
  "with RBF gamma = %0.4f " % gamma)

print("\n============================ ")
print("\n1. Using FFT - tricky RBF ")
rnd = np.random.RandomState(0)
prune_mask = prune_FFT_1(X, n_to_keep, rnd, gamma)
print("\nPrune mask: ")
print(prune_mask)

print("\n============================ ")
print("\n2. Using FFT - clear RBF ")
rnd = np.random.RandomState(0)
prune_mask = prune_FFT_2(X, n_to_keep, rnd, gamma)
print("\nPrune mask: ")
print(prune_mask)

print("\n============================ ")
print("\n3. Using KRA - full memory K ")
prune_mask = prune_KRA_1(X, n_to_keep, gamma)
print("\nPrune mask: ")
print(prune_mask)

print("\n============================ ")
print("\n4. Using KRA - low memory K ")
prune_mask = prune_KRA_2(X, n_to_keep, gamma)
print("\nPrune mask: ")
print(prune_mask)
print("\n============================ ")

print("\nEnd demo ")

An Extra Trees Regressor (“Extremely Randomized Trees”) is an ensemble machine learning algorithm that builds multiple randomized decision trees to predict continuous numerical values. It differs from Random Forest regressions by introducing greater randomization during node splitting and utilizing the entire original dataset to train each tree rather than using samples of the original dataset.

One weekend I figured I’d implement Extra Trees Regression, from scratch, using Python and NumPy. The process took me a few hours, but was manageable because I have quite a bit of experience with tree-based systems.

The output of my demo:

Extra trees regression scratch Python NumPy

Loading synthetic train (200), test (40) data
Done

First three X predictors:
[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]
[ 0.0776 -0.1616  0.3704 -0.5911  0.7562]
[-0.9452  0.3409 -0.1654  0.1174 -0.7192]

First three y targets:
0.4840
0.1568
0.8054

Setting n_estimators = 200
Setting max_depth = 6
Setting max_features = 4

Accuracy on train (within 0.10) = 0.7700
Accuracy on test (within 0.10) = 0.7250

MSE on train = 0.0006
MSE on test = 0.0012

Predicting for:
[[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]]
Predicted y = 0.4930

SCIKIT:

Accuracy on train (within 0.10) = 0.7750
Accuracy on test (within 0.10) = 0.6750

MSE on train = 0.0006
MSE on test = 0.0012

Done

The data is synthetic. It was generated by a neural network with random weights and biases. There are five predictor values and one value to predict. There are 200 training items and 40 test items.

Extra Trees Regression is an ensemble technique. For my demo, I used a collection of 200 base decision trees. The final prediction is the average of the 200 base trees. The base trees are exactly like standard decision trees, except, during construction, at each node, a base Extra Tree picks one value at random from each active column, and then uses the one value that gives the best split — a separation of rows of data that has the lowest variance (so that the left child nodes and right child nodes have relatively similar target y values).

The ideas are suprisingly subtle and implementation is extremely tricky.

The output of my from-scratch demo is very close to the output from the scikit-learn ExtraTreesRegressor module. The output wasn’t expected to be identical because Extra Trees Regression is fundamentally random.

Extra Trees Regression has no super powers — it’s just math and computer science. When I was a young man, I loved, and learned to read from, comic books of the early 1960s. I especially loved heros who had no special powers — they relied on just their intelligence — Batman, Green Arrow, Adam Strange, and Iron Man.

Right: In issue #48 (December 1963), the suit changed to red and gold. Iron Man got his own title in May 1968.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. (My blog editor chokes on symbols).

# extra_trees_regression_scratch.py

import numpy as np

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed', linewidth=120)

# ===========================================================

class MyExtraTreesRegressor:
  def __init__(self, n_estimators=100, max_depth=3,
    min_samples_split=2, max_features='auto',
    RandomState=0):

    self.n_estimators = n_estimators
    self.max_depth = max_depth
    self.min_samples_split = min_samples_split
    self.max_features = max_features
    self.RandomState = RandomState
    self.trees = []
    self.rng = np.random.RandomState(RandomState)

  def fit(self, train_X, train_y):
    X = np.asarray(train_X)
    n_samples, n_features = X.shape

    # map scikit-learn standard max_features to your
    # n_split_cols parameter
    if self.max_features == 'auto' or self.max_features is None:
      n_split_cols = n_features
    elif self.max_features == 'sqrt':
      n_split_cols = int(np.sqrt(n_features))
    elif self.max_features == 'log2':
      n_split_cols = int(np.log2(n_features))
    else:
      n_split_cols = int(self.max_features)
    n_split_cols = max(1, min(n_split_cols, n_features))

    self.trees = []
    for _ in range(self.n_estimators):
      tree_seed = self.rng.randint(0, 1000000)
      
      # instantiate our specialized base tree 
      tree = MyExtraTreesBaseTree(
        max_depth=self.max_depth,
        min_samples=self.min_samples_split,
        min_leaf=1,
        n_split_cols=n_split_cols,
        save_rows=False,
        seed=tree_seed
      )
      tree.fit(train_X, train_y)
      self.trees.append(tree)
    return self

  def predict(self, X):
    X = np.asarray(X)
    preds = np.array([tree.predict(X) for tree in self.trees])
    return np.mean(preds, axis=0)


# ===========================================================

# -----------------------------------------------------------

def accuracy(model, data_X, data_y, pct_close):
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    y_pred = model.predict(x)[0]

    if np.abs(y - y_pred) "lt" np.abs(y * pct_close):
      n_correct += 1
    else: 
      n_wrong += 1
  # print("Correct = " + str(n_correct))
  # print("Wrong   = " + str(n_wrong))
  return n_correct / (n_correct + n_wrong)

# -----------------------------------------------------------

def MSE(model, data_X, data_y):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    y_pred = model.predict(x)[0]
    # print(y_pred); input()
    sum += (y - y_pred) * (y - y_pred)

  return sum / n

# -----------------------------------------------------------

# ===========================================================

class MyExtraTreesBaseTree:

  # identical to MyDecisionTreeRegressor except best_split()

  def __init__(self, max_depth=3, min_samples=2, 
    min_leaf=1, n_split_cols=-1, save_rows=False, 
    seed=0):

    self.max_depth = max_depth
    self.min_samples = min_samples   # to split
    self.min_leaf = min_leaf         # after split
    self.n_split_cols = n_split_cols # -1 all cols
    self.save_rows = save_rows
    self.rnd = np.random.RandomState(seed)

    self.tree = []  # list: no recursive, no ptrs
    num_nodes = 2**(max_depth + 1) - 1
    for i in range(num_nodes):
      self.tree.append(None)

    self.train_X = None
    self.train_y = None

  # ...............................................

  class Node:
    def __init__(self):
      self.id = -1
      self.col_idx = -1  # aka split col
      self.thresh = 0.0  # aka split val
      self.left = -1
      self.right = -1
      self.value = 0.0   # predicted y
      self.is_leaf = False
      self.rows = []     # associated rows in train

  # ...............................................

  def fit(self, train_X, train_y):
    # aka train(), build_tree()
    # self.train_X = np.array(train_X)  # safety
    # self.train_y = np.array(train_y)
    self.train_X = train_X
    self.train_y = train_y
    
    max_id = 2**(self.max_depth+1) - 2  # md = 3, m_id = 14
    max_start_id = 2**self.max_depth - 1  # md = 3, ms = 7

    # prep root node
    all_rows = []
    for i in range(len(self.train_X)):  # or use arange
      all_rows.append(i)
    grand_mean = self.tree_target_mean(all_rows)
    
    root = self.Node()
    root.id = 0
    root.value = grand_mean
    root.is_leaf = False
    root.rows = all_rows
    self.tree[0] = root

    for i in range(len(self.tree)):  # each node
      curr_node = self.tree[i]       # convenience
    
      if curr_node is None: continue

      # if node too deep to have children
      #  or not enough rows to split, leave children alone
      if curr_node.id "gte" max_start_id or \
        len(curr_node.rows) "lt" self.min_samples:
        curr_node.is_leaf = True
        continue

      # try to split curr node
      col_idx, split_val = self.best_split(curr_node.rows)
      if col_idx == -1:  # bad split
        curr_node.is_leaf = True
        curr_node.left = -1
        curr_node.right = -1
        continue
  
      # got good split info
      curr_node.col_idx = col_idx
      curr_node.thresh = split_val

      # make rows for the children
      left_idxs = []
      right_idxs = []
      for k in range(len(curr_node.rows)):
        r = curr_node.rows[k]
        if self.train_X[r][col_idx] "lte" split_val:
          left_idxs.append(r)
        else:
          right_idxs.append(r)

      # make left child
      left_id = curr_node.id * 2 + 1
      if left_id "lte" max_id and len(left_idxs) "gte" \
        self.min_leaf:
        curr_node.left = left_id
        left_node = self.Node()
        left_node.id = left_id
        left_node.rows = left_idxs
        left_node.value = \
          self.tree_target_mean(left_node.rows)
        self.tree[left_id] = left_node
      else:
        curr_node.left = -1

      # make right child
      right_id = curr_node.id * 2 + 2
      if right_id "lte" max_id and len(right_idxs) "gte" \
        self.min_leaf:
        curr_node.right = right_id
        right_node = self.Node()
        right_node.id = right_id
        right_node.rows = right_idxs
        right_node.value = \
          self.tree_target_mean(right_node.rows)
        self.tree[right_id] = right_node
      else:
        curr_node.right = -1

      if curr_node.left == -1 and curr_node.right == -1:
        curr_node.is_leaf = True

    # zap away rows to save space if in ensemble
    if self.save_rows == False:
      for k in range(len(self.tree)):
        if self.tree[k] is not None:
          self.tree[k].rows = None

  # ---------------------------------------------------------

  def predict_one(self, x):
    # x is a vector
    curr_idx = 0
    last_valid_value = 0.0
        
    while curr_idx != -1 and curr_idx "lt" len(self.tree):
      curr_node = self.tree[curr_idx]
      
      if curr_node is None: break  # safety check
                
      last_valid_value = curr_node.value
      if curr_node.is_leaf == True: break

      if x[curr_node.col_idx] "lte" curr_node.thresh:
        curr_idx = curr_node.left
      else:
        curr_idx = curr_node.right

    return last_valid_value

  # ---------------------------------------------------------

  def predict(self, X):
    # X is a matrix
    n = len(X)
    result = np.zeros(n)
    for i in range(n):
      result[i] = self.predict_one(X[i])
    return result

  # ---------------------------------------------------------

  def best_split(self, rows):  # Extra Trees version
    best_col_idx = -1
    best_thresh = 0.0
    best_var = float('inf')
    n_rows = len(rows)
    n_cols = len(self.train_X[0])

    # Fisher-Yates shuffle to get candidate features
    col_indices = np.arange(n_cols)
    for i in range(n_cols - 1):
      ri = self.rnd.randint(i, n_cols)
      tmp = col_indices[i]
      col_indices[i] = col_indices[ri]
      col_indices[ri] = tmp

    n_cols_to_use = min(self.n_split_cols, n_cols) \
      if self.n_split_cols != -1 else n_cols
    active_cols = col_indices[0:n_cols_to_use]

    # evaluate a threshold for each active column
    for col_idx in active_cols:
      feature_vals = self.train_X[rows, col_idx]
      min_v = np.min(feature_vals)
      max_v = np.max(feature_vals)
      
      # if all feature values are identical, 
      # this column cannot be split
      if min_v == max_v:
        continue
      
      # EXTRA TREES MECHANIC: 
      # select a single threshold completely at random
      thresh = self.rnd.uniform(min_v, max_v)
      
      left_count = 0
      right_count = 0
      left_sum_y = 0.0
      left_sum_sq_y = 0.0
      
      # count and sum target values to compute
      # variance reduction
      for r in rows:
        y_curr = self.train_y[r]
        if self.train_X[r, col_idx] "lte" thresh:
          left_count += 1
          left_sum_y += y_curr
          left_sum_sq_y += y_curr * y_curr
        else:
          right_count += 1

      # enforce min_leaf boundaries
      if left_count "lt" self.min_leaf or \
        right_count "lt" self.min_leaf:
        continue

      # isolate right side aggregates dynamically
      # using total sums
      total_sum_y = np.sum(self.train_y[rows])
      total_sum_sq_y = np.sum(self.train_y[rows] ** 2)
      
      right_sum_y = total_sum_y - left_sum_y
      right_sum_sq_y = total_sum_sq_y - left_sum_sq_y

      # variance shortcut formula: E[X^2] - (E[X])^2
      left_var = (left_sum_sq_y / left_count) - \
        ((left_sum_y / left_count) ** 2)
      right_var = (right_sum_sq_y / right_count) - \
        ((right_sum_y / right_count) ** 2)
      
      if left_var "lt" 0.0: left_var = 0.0
      if right_var "lt" 0.0: right_var = 0.0
      
      weighted_var = ((left_count * left_var) + \
        (right_count * right_var)) / n_rows
      
      if weighted_var "lt" best_var:
        best_var = weighted_var
        best_col_idx = col_idx
        best_thresh = thresh

    return best_col_idx, best_thresh

  # ---------------------------------------------------------

  def tree_target_mean(self, rows):
    if rows is None or len(rows) == 0: return 0.0
    sum = 0.0
    for i in range(len(rows)):
      r = rows[i]
      sum += self.train_y[r]
    return sum / len(rows)

  # ---------------------------------------------------------

# ===========================================================
# ===========================================================

print("\nExtra trees regression scratch Python NumPy ")

print("\nLoading synthetic train (200), test (40) data ")
train_file = ".\\Data\\synthetic_train_200.txt"
train_X = np.loadtxt(train_file, comments="#",
  usecols=[0,1,2,3,4],
  delimiter=",",  dtype=np.float64)
train_y = np.loadtxt(train_file, comments="#", usecols=5,
  delimiter=",",  dtype=np.float64)

test_file = ".\\Data\\synthetic_test_40.txt"
test_X = np.loadtxt(test_file, comments="#",
  usecols=[0,1,2,3,4],
  delimiter=",",  dtype=np.float64)
test_y = np.loadtxt(test_file, comments="#", usecols=5,
  delimiter=",",  dtype=np.float64)
print("Done ")

print("\nFirst three X predictors: ")
for i in range(3):
  print(train_X[i])
print("\nFirst three y targets: ")
for i in range(3):
  print("%0.4f" % train_y[i])

n_ests = 200
max_d = 6
max_f = 4

print("\nSetting n_estimators = " + str(n_ests))
print("Setting max_depth = " + str(max_d))
print("Setting max_features = " + str(max_f))

model = MyExtraTreesRegressor(n_estimators=n_ests,
 max_depth=max_d, min_samples_split=2, 
 max_features=max_f, RandomState=0)

model.fit(train_X, train_y)

acc_train = accuracy(model, train_X, train_y, 0.10)
print("\nAccuracy on train (within 0.10) = \
%0.4f " % acc_train)
acc_test = accuracy(model, test_X, test_y, 0.10)
print("Accuracy on test (within 0.10) = \
%0.4f " % acc_test)

mse_train = MSE(model, train_X, train_y)
print("\nMSE on train = %0.4f " % mse_train)
mse_test = MSE(model, test_X, test_y)
print("MSE on test = %0.4f " % mse_test)

x = train_X[0].reshape(1,-1)
print("\nPredicting for: ")
print(x)
y_pred = model.predict(x)[0]
print("Predicted y = %0.4f " % y_pred)

# ==== SCIKIT ===========

from sklearn.ensemble import ExtraTreesRegressor

print("\nSCIKIT: ")

model = ExtraTreesRegressor(n_estimators=n_ests, 
  max_depth=max_d, max_features=max_f, random_state=0)
model.fit(train_X, train_y)

acc_train = accuracy(model, train_X, train_y, 0.10)
print("\nAccuracy on train (within 0.10) = \
%0.4f " % acc_train)
acc_test = accuracy(model, test_X, test_y, 0.10)
print("Accuracy on test (within 0.10) = \
%0.4f " % acc_test)

mse_train = MSE(model, train_X, train_y)
print("\nMSE on train = %0.4f " % mse_train)
mse_test = MSE(model, test_X, test_y)
print("MSE on test = %0.4f " % mse_test)

print("\nDone ")

Training data:

# synthetic_train_200.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996
-0.9434, -0.5076,  0.7201,  0.0777,  0.1056,  0.5664
 0.9392,  0.1221, -0.9627,  0.6013, -0.5341,  0.1533
 0.6142, -0.2243,  0.7271,  0.4942,  0.1125,  0.1661
 0.4260,  0.1194, -0.9749, -0.8561,  0.9346,  0.2230
 0.1362, -0.5934, -0.4953,  0.4877, -0.6091,  0.3810
 0.6937, -0.5203, -0.0125,  0.2399,  0.6580,  0.1460
-0.6864, -0.9628, -0.8600, -0.0273,  0.2127,  0.5387
 0.9772,  0.1595, -0.2397,  0.1019,  0.4907,  0.1611
 0.3385, -0.4702, -0.8673, -0.2598,  0.2594,  0.2270
-0.8669, -0.4794,  0.6095, -0.6131,  0.2789,  0.4700
 0.0493,  0.8496, -0.4734, -0.8681,  0.4701,  0.3516
 0.8639, -0.9721, -0.5313,  0.2336,  0.8980,  0.1412
 0.9004,  0.1133,  0.8312,  0.2831, -0.2200,  0.1782
 0.0991,  0.8524,  0.8375, -0.2102,  0.9265,  0.2150
-0.6521, -0.7473, -0.7298,  0.0113, -0.9570,  0.7422
 0.6190, -0.3105,  0.8802,  0.1640,  0.7577,  0.1056
 0.6895,  0.8108, -0.0802,  0.0927,  0.5972,  0.2214
 0.1982, -0.9689,  0.1870, -0.1326,  0.6147,  0.1310
-0.3695,  0.7858,  0.1557, -0.6320,  0.5759,  0.3773
-0.1596,  0.3581,  0.8372, -0.9992,  0.9535,  0.2071
-0.2468,  0.9476,  0.2094,  0.6577,  0.1494,  0.4132
 0.1737,  0.5000,  0.7166,  0.5102,  0.3961,  0.2611
 0.7290, -0.3546,  0.3416, -0.0983, -0.2358,  0.1332
-0.3652,  0.2438, -0.1395,  0.9476,  0.3556,  0.4170
-0.6029, -0.1466, -0.3133,  0.5953,  0.7600,  0.4334
-0.4596, -0.4953,  0.7098,  0.0554,  0.6043,  0.2775
 0.1450,  0.4663,  0.0380,  0.5418,  0.1377,  0.2931
-0.8636, -0.2442, -0.8407,  0.9656, -0.6368,  0.7429
 0.6237,  0.7499,  0.3768,  0.1390, -0.6781,  0.2185
-0.5499,  0.1850, -0.3755,  0.8326,  0.8193,  0.4399
-0.4858, -0.7782, -0.6141, -0.0008,  0.4572,  0.4197
 0.7033, -0.1683,  0.2334, -0.5327, -0.7961,  0.1776
 0.0317, -0.0457, -0.6947,  0.2436,  0.0880,  0.3345
 0.5031, -0.5559,  0.0387,  0.5706, -0.9553,  0.3107
-0.3513,  0.7458,  0.6894,  0.0769,  0.7332,  0.3170
 0.2205,  0.5992, -0.9309,  0.5405,  0.4635,  0.3532
-0.4806, -0.4859,  0.2646, -0.3094,  0.5932,  0.3202
 0.9809, -0.3995, -0.7140,  0.8026,  0.0831,  0.1600
 0.9495,  0.2732,  0.9878,  0.0921,  0.0529,  0.1289
-0.9476, -0.6792,  0.4913, -0.9392, -0.2669,  0.5966
 0.7247,  0.3854,  0.3819, -0.6227, -0.1162,  0.1550
-0.5922, -0.5045, -0.4757,  0.5003, -0.0860,  0.5863
-0.8861,  0.0170, -0.5761,  0.5972, -0.4053,  0.7301
 0.6877, -0.2380,  0.4997,  0.0223,  0.0819,  0.1404
 0.9189,  0.6079, -0.9354,  0.4188, -0.0700,  0.1907
-0.1428, -0.7820,  0.2676,  0.6059,  0.3936,  0.2790
 0.5324, -0.3151,  0.6917, -0.1425,  0.6480,  0.1071
-0.8432, -0.9633, -0.8666, -0.0828, -0.7733,  0.7784
-0.9444,  0.5097, -0.2103,  0.4939, -0.0952,  0.6787
-0.0520,  0.6063, -0.1952,  0.8094, -0.9259,  0.4836
 0.5477, -0.7487,  0.2370, -0.9793,  0.0773,  0.1241
 0.2450,  0.8116,  0.9799,  0.4222,  0.4636,  0.2355
 0.8186, -0.1983, -0.5003, -0.6531, -0.7611,  0.1511
-0.4714,  0.6382, -0.3788,  0.9648, -0.4667,  0.5950
 0.0673, -0.3711,  0.8215, -0.2669, -0.1328,  0.2677
-0.9381,  0.4338,  0.7820, -0.9454,  0.0441,  0.5518
-0.3480,  0.7190,  0.1170,  0.3805, -0.0943,  0.4724
-0.9813,  0.1535, -0.3771,  0.0345,  0.8328,  0.5438
-0.1471, -0.5052, -0.2574,  0.8637,  0.8737,  0.3042
-0.5454, -0.3712, -0.6505,  0.2142, -0.1728,  0.5783
 0.6327, -0.6297,  0.4038, -0.5193,  0.1484,  0.1153
-0.5424,  0.3282, -0.0055,  0.0380, -0.6506,  0.6613
 0.1414,  0.9935,  0.6337,  0.1887,  0.9520,  0.2540
-0.9351, -0.8128, -0.8693, -0.0965, -0.2491,  0.7353
 0.9507, -0.6640,  0.9456,  0.5349,  0.6485,  0.1059
-0.0462, -0.9737, -0.2940, -0.0159,  0.4602,  0.2606
-0.0627, -0.0852, -0.7247, -0.9782,  0.5166,  0.2977
 0.0478,  0.5098, -0.0723, -0.7504, -0.3750,  0.3335
 0.0090,  0.3477,  0.5403, -0.7393, -0.9542,  0.4415
-0.9748,  0.3449,  0.3736, -0.1015,  0.8296,  0.4358
 0.2887, -0.9895, -0.0311,  0.7186,  0.6608,  0.2057
 0.1570, -0.4518,  0.1211,  0.3435, -0.2951,  0.3244
 0.7117, -0.6099,  0.4946, -0.4208,  0.5476,  0.1096
-0.2929, -0.5726,  0.5346, -0.3827,  0.4665,  0.2465
 0.4889, -0.5572, -0.5718, -0.6021, -0.7150,  0.2163
-0.7782,  0.3491,  0.5996, -0.8389, -0.5366,  0.6516
-0.5847,  0.8347,  0.4226,  0.1078, -0.3910,  0.6134
 0.8469,  0.4121, -0.0439, -0.7476,  0.9521,  0.1571
-0.6803, -0.5948, -0.1376, -0.1916, -0.7065,  0.7156
 0.2878,  0.5086, -0.5785,  0.2019,  0.4979,  0.2980
 0.2764,  0.1943, -0.4090,  0.4632,  0.8906,  0.2960
-0.8877,  0.6705, -0.6155, -0.2098, -0.3998,  0.7107
-0.8398,  0.8093, -0.2597,  0.0614, -0.0118,  0.6502
-0.8476,  0.0158, -0.4769, -0.2859, -0.7839,  0.7715
 0.5751, -0.7868,  0.9714, -0.6457,  0.1448,  0.1175
 0.4802, -0.7001,  0.1022, -0.5668,  0.5184,  0.1090
 0.4458, -0.6469,  0.7239, -0.9604,  0.7205,  0.0779
 0.5175,  0.4339,  0.9747, -0.4438, -0.9924,  0.2879
 0.8678,  0.7158,  0.4577,  0.0334,  0.4139,  0.1678
 0.5406,  0.5012,  0.2264, -0.1963,  0.3946,  0.2088
-0.9938,  0.5498,  0.7928, -0.5214, -0.7585,  0.7687
 0.7661,  0.0863, -0.4266, -0.7233, -0.4197,  0.1466
 0.2277, -0.3517, -0.0853, -0.1118,  0.6563,  0.1767
 0.3499, -0.5570, -0.0655, -0.3705,  0.2537,  0.1632
 0.7547, -0.1046,  0.5689, -0.0861,  0.3125,  0.1257
 0.8186,  0.2110,  0.5335,  0.0094, -0.0039,  0.1391
 0.6858, -0.8644,  0.1465,  0.8855,  0.0357,  0.1845
-0.4967,  0.4015,  0.0805,  0.8977,  0.2487,  0.4663
 0.6760, -0.9841,  0.9787, -0.8446, -0.3557,  0.1509
-0.1203, -0.4885,  0.6054, -0.0443, -0.7313,  0.4854
 0.8557,  0.7919, -0.0169,  0.7134, -0.1628,  0.2002
 0.0115, -0.6209,  0.9300, -0.4116, -0.7931,  0.4052
-0.7114, -0.9718,  0.4319,  0.1290,  0.5892,  0.3661
 0.3915,  0.5557, -0.1870,  0.2955, -0.6404,  0.2954
-0.3564, -0.6548, -0.1827, -0.5172, -0.1862,  0.4622
 0.2392, -0.4959,  0.5857, -0.1341, -0.2850,  0.2470
-0.3394,  0.3947, -0.4627,  0.6166, -0.4094,  0.5325
 0.7107,  0.7768, -0.6312,  0.1707,  0.7964,  0.2757
-0.1078,  0.8437, -0.4420,  0.2177,  0.3649,  0.4028
-0.3139,  0.5595, -0.6505, -0.3161, -0.7108,  0.5546
 0.4335,  0.3986,  0.3770, -0.4932,  0.3847,  0.1810
-0.2562, -0.2894, -0.8847,  0.2633,  0.4146,  0.4036
 0.2272,  0.2966, -0.6601, -0.7011,  0.0284,  0.2778
-0.0743, -0.1421, -0.0054, -0.6770, -0.3151,  0.3597
-0.4762,  0.6891,  0.6007, -0.1467,  0.2140,  0.4266
-0.4061,  0.7193,  0.3432,  0.2669, -0.7505,  0.6147
-0.0588,  0.9731,  0.8966,  0.2902, -0.6966,  0.4955
-0.0627, -0.1439,  0.1985,  0.6999,  0.5022,  0.3077
 0.1587,  0.8494, -0.8705,  0.9827, -0.8940,  0.4263
-0.7850,  0.2473, -0.9040, -0.4308, -0.8779,  0.7199
 0.4070,  0.3369, -0.2428, -0.6236,  0.4940,  0.2215
-0.0242,  0.0513, -0.9430,  0.2885, -0.2987,  0.3947
-0.5416, -0.1322, -0.2351, -0.0604,  0.9590,  0.3683
 0.1055,  0.7783, -0.2901, -0.5090,  0.8220,  0.2984
-0.9129,  0.9015,  0.1128, -0.2473,  0.9901,  0.4776
-0.9378,  0.1424, -0.6391,  0.2619,  0.9618,  0.5368
 0.7498, -0.0963,  0.4169,  0.5549, -0.0103,  0.1614
-0.2612, -0.7156,  0.4538, -0.0460, -0.1022,  0.3717
 0.7720,  0.0552, -0.1818, -0.4622, -0.8560,  0.1685
-0.4177,  0.0070,  0.9319, -0.7812,  0.3461,  0.3052
-0.0001,  0.5542, -0.7128, -0.8336, -0.2016,  0.3803
 0.5356, -0.4194, -0.5662, -0.9666, -0.2027,  0.1776
-0.2378,  0.3187, -0.8582, -0.6948, -0.9668,  0.5474
-0.1947, -0.3579,  0.1158,  0.9869,  0.6690,  0.2992
 0.3992,  0.8365, -0.9205, -0.8593, -0.0520,  0.3154
-0.0209,  0.0793,  0.7905, -0.1067,  0.7541,  0.1864
-0.4928, -0.4524, -0.3433,  0.0951, -0.5597,  0.6261
-0.8118,  0.7404, -0.5263, -0.2280,  0.1431,  0.6349
 0.0516, -0.8480,  0.7483,  0.9023,  0.6250,  0.1959
-0.3212,  0.1093,  0.9488, -0.3766,  0.3376,  0.2735
-0.3481,  0.5490, -0.3484,  0.7797,  0.5034,  0.4379
-0.5785, -0.9170, -0.3563, -0.9258,  0.3877,  0.4121
 0.3407, -0.1391,  0.5356,  0.0720, -0.9203,  0.3458
-0.3287, -0.8954,  0.2102,  0.0241,  0.2349,  0.3247
-0.1353,  0.6954, -0.0919, -0.9692,  0.7461,  0.3338
 0.9036, -0.8982, -0.5299, -0.8733, -0.1567,  0.1187
 0.7277, -0.8368, -0.0538, -0.7489,  0.5458,  0.0830
 0.9049,  0.8878,  0.2279,  0.9470, -0.3103,  0.2194
 0.7957, -0.1308, -0.5284,  0.8817,  0.3684,  0.2172
 0.4647, -0.4931,  0.2010,  0.6292, -0.8918,  0.3371
-0.7390,  0.6849,  0.2367,  0.0626, -0.5034,  0.7039
-0.1567, -0.8711,  0.7940, -0.5932,  0.6525,  0.1710
 0.7635, -0.0265,  0.1969,  0.0545,  0.2496,  0.1445
 0.7675,  0.1354, -0.7698, -0.5460,  0.1920,  0.1728
-0.5211, -0.7372, -0.6763,  0.6897,  0.2044,  0.5217
 0.1913,  0.1980,  0.2314, -0.8816,  0.5006,  0.1998
 0.8964,  0.0694, -0.6149,  0.5059, -0.9854,  0.1825
 0.1767,  0.7104,  0.2093,  0.6452,  0.7590,  0.2832
-0.3580, -0.7541,  0.4426, -0.1193, -0.7465,  0.5657
-0.5996,  0.5766, -0.9758, -0.3933, -0.9572,  0.6800
 0.9950,  0.1641, -0.4132,  0.8579,  0.0142,  0.2003
-0.4717, -0.3894, -0.2567, -0.5111,  0.1691,  0.4266
 0.3917, -0.8561,  0.9422,  0.5061,  0.6123,  0.1212
-0.0366, -0.1087,  0.3449, -0.1025,  0.4086,  0.2475
 0.3633,  0.3943,  0.2372, -0.6980,  0.5216,  0.1925
-0.5325, -0.6466, -0.2178, -0.3589,  0.6310,  0.3568
 0.2271,  0.5200, -0.1447, -0.8011, -0.7699,  0.3128
 0.6415,  0.1993,  0.3777, -0.0178, -0.8237,  0.2181
-0.5298, -0.0768, -0.6028, -0.9490,  0.4588,  0.4356
 0.6870, -0.1431,  0.7294,  0.3141,  0.1621,  0.1632
-0.5985,  0.0591,  0.7889, -0.3900,  0.7419,  0.2945
 0.3661,  0.7984, -0.8486,  0.7572, -0.6183,  0.3449
 0.6995,  0.3342, -0.3113, -0.6972,  0.2707,  0.1712
 0.2565,  0.9126,  0.1798, -0.6043, -0.1413,  0.2893
-0.3265,  0.9839, -0.2395,  0.9854,  0.0376,  0.4770
 0.2690, -0.1722,  0.9818,  0.8599, -0.7015,  0.3954
-0.2102, -0.0768,  0.1219,  0.5607, -0.0256,  0.3949
 0.8216, -0.9555,  0.6422, -0.6231,  0.3715,  0.0801
-0.2896,  0.9484, -0.7545, -0.6249,  0.7789,  0.4370
-0.9985, -0.5448, -0.7092, -0.5931,  0.7926,  0.5402

Test data:

# synthetic_test_40.txt
#
 0.7462,  0.4006, -0.0590,  0.6543, -0.0083,  0.1935
 0.8495, -0.2260, -0.0142, -0.4911,  0.7699,  0.1078
-0.2335, -0.4049,  0.4352, -0.6183, -0.7636,  0.5088
 0.1810, -0.5142,  0.2465,  0.2767, -0.3449,  0.3136
-0.8650,  0.7611, -0.0801,  0.5277, -0.4922,  0.7140
-0.2358, -0.7466, -0.5115, -0.8413, -0.3943,  0.4533
 0.4834,  0.2300,  0.3448, -0.9832,  0.3568,  0.1360
-0.6502, -0.6300,  0.6885,  0.9652,  0.8275,  0.3046
-0.3053,  0.5604,  0.0929,  0.6329, -0.0325,  0.4756
-0.7995,  0.0740, -0.2680,  0.2086,  0.9176,  0.4565
-0.2144, -0.2141,  0.5813,  0.2902, -0.2122,  0.4119
-0.7278, -0.0987, -0.3312, -0.5641,  0.8515,  0.4438
 0.3793,  0.1976,  0.4933,  0.0839,  0.4011,  0.1905
-0.8568,  0.9573, -0.5272,  0.3212, -0.8207,  0.7415
-0.5785,  0.0056, -0.7901, -0.2223,  0.0760,  0.5551
 0.0735, -0.2188,  0.3925,  0.3570,  0.3746,  0.2191
 0.1230, -0.2838,  0.2262,  0.8715,  0.1938,  0.2878
 0.4792, -0.9248,  0.5295,  0.0366, -0.9894,  0.3149
-0.4456,  0.0697,  0.5359, -0.8938,  0.0981,  0.3879
 0.8629, -0.8505, -0.4464,  0.8385,  0.5300,  0.1769
 0.1995,  0.6659,  0.7921,  0.9454,  0.9970,  0.2330
-0.0249, -0.3066, -0.2927, -0.4923,  0.8220,  0.2437
 0.4513, -0.9481, -0.0770, -0.4374, -0.9421,  0.2879
-0.3405,  0.5931, -0.3507, -0.3842,  0.8562,  0.3987
 0.9538,  0.0471,  0.9039,  0.7760,  0.0361,  0.1706
-0.0887,  0.2104,  0.9808,  0.5478, -0.3314,  0.4128
-0.8220, -0.6302,  0.0537, -0.1658,  0.6013,  0.4306
-0.4123, -0.2880,  0.9074, -0.0461, -0.4435,  0.5144
 0.0060,  0.2867, -0.7775,  0.5161,  0.7039,  0.3599
-0.7968, -0.5484,  0.9426, -0.4308,  0.8148,  0.2979
 0.7811,  0.8450, -0.6877,  0.7594,  0.2640,  0.2362
-0.6802, -0.1113, -0.8325, -0.6694, -0.6056,  0.6544
 0.3821,  0.1476,  0.7466, -0.5107,  0.2592,  0.1648
 0.7265,  0.9683, -0.9803, -0.4943, -0.5523,  0.2454
-0.9049, -0.9797, -0.0196, -0.9090, -0.4433,  0.6447
-0.4607,  0.1811, -0.2389,  0.4050, -0.0078,  0.5229
 0.2664, -0.2932, -0.4259, -0.7336,  0.8742,  0.1834
-0.4507,  0.1029, -0.6294, -0.1158, -0.6294,  0.6081
 0.8948, -0.0124,  0.9278,  0.2899, -0.0314,  0.1534
-0.1323, -0.8813, -0.0146, -0.0697,  0.6135,  0.2386

In machine learning, if training data is multicollinear, the interpretability of the resulting model will likely be poor (the details are tricky and out of the scope of this post). Multicollinear data has two or more rows that are highly correlated mathematically. For example, if a set of training data has a column (predictor variable) house-size (in square feet), and another column house-price, those two columns are likely highly correlated with each other.

Note: Multicollinear training data is usually not a problem for neural network regression, and not a problem for tree-based regression (random forest, gradient boosting).

The most common way to analyze training data for multicollinearity is to compute the VIF (variance inflation factor) for each column of the data.

VIF is a value between 1.0 and positive infinity (in weird scenarios, a VIF value could be less than 1.0). Briefly, if all column VIF values are less than about 7.0, the data is probably OK in terms of multicollinearity.

if VIF is close to 1.0, the column is not correlated with other cols
if VIF between 1.0 and 5.0, column is mildly correlated
if VIF between 5.0 and 10.0, column is highly correlated
if VIF greater than 10.0, column is extremely correlated

To compute the VIF for a specified column of training data, you use the specified column as the dependent y variable, and use the remaining columns as the independent predictor variables, and compute a linear regression model, and then compute the R2 (coefficient of determination) for the model. The VIF value for the column is 1.0 / (1.0 – R2).

Suppose that you have a set of training data X predictor values, and you use some column c as the dependent y variable, and all the other columns as predictors for c. After training a linear regression model, you compute R2 and it is 0.90 — which means column c is predicted very well by the other columns. The VIF value for column c is 1.0 / (1.0 – R2) = 1.0 / 0.10 = 10.0 which is large which is bad because column c is a linear combination of the other columns — the data is somewhat multicollinear.

Now, with the same setup, suppose R2 is 0.20 — which means column c cannot be predicted well by the other columns. The VIF value is 1.0 / (1.0 – 0.20) = 1.0 / 0.8 = 1.25 which is a small value, which is good, because column c is not a linear combination of the other columns, and therefore the data is not multicollinear.

Some time ago, I put together a demo using Python NumPy and the scikit library. This was easy because NumPy and scikit had all the components available.

I decided to entertain myself by refactoring my Python demo to the C# language. This required a bit of effort because I had to implement a no-touch linear regression class, plus several helper functions, such as a program-defined MatWithoutColumn() to mimic the built-in numpy np.delete() function.

I created two datasets. The first data set has five columns of predictors, followed by a column of target y values. The data is “normal” in the sense that there’s no multicollinearity. There are 20 items. It looks like:

-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
. . .

The second dataset is highly multicollinear, where the third column is 2 times the first column, plus the second column, plus a small random value between 0.000 and 0.001. It looks like:

-0.1660,  0.4406,  0.1096, -0.3953, -0.7065, 0.4840
 0.0776, -0.1616, -0.0045, -0.5911,  0.7562, 0.1568
-0.9452,  0.3409, -1.5482,  0.1174, -0.7192, 0.8054
. . .

The output of my C# VIF demo program is:

Begin variance inflation factor (VIF) demo using C#

Loading synthetic (20) normal data

First three lines:
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065
   0.0776  -0.1616   0.3704  -0.5911   0.7562
  -0.9452   0.3409  -0.1654   0.1174  -0.7192

Begin VIF analysis
col =  0  |  vif = 1.1979
col =  1  |  vif = 1.4590
col =  2  |  vif = 1.2345
col =  3  |  vif = 1.3025
col =  4  |  vif = 1.2120

Loading synthetic (20) highly multicollinear data
(col[2] = 2.0 * col[0] + col[1] + noise)

First three lines:
  -0.1660   0.4406   0.1096  -0.3953  -0.7065
   0.0776  -0.1616  -0.0045  -0.5911   0.7562
  -0.9452   0.3409  -1.5482   0.1174  -0.7192

Begin VIF analysis
col =  0  |  vif = 25101680.1917
col =  1  |  vif = 5769710.7907
col =  2  |  vif = 30816932.6001
col =  3  |  vif = 1.2937
col =  4  |  vif = 1.1879

End demo

As expected, the first dataset didn’t have any bad VIF values. But the VIF values for the second dataset show that columns [0], [1], [2] are highly correlated.

For my internal linear regressor, I used SGD training with some default parameters and auto-early-exit, so that the regressor just works without tuning. While experimenting, I noticed that SGD training has a lot of trouble dealing with multicollinear data — which is the whole point of checking for multicollinear data in the first place.

An interesting exploration.

In machine learning, you don’t want a relationship between two columns in your training data. But in science fiction movies, you absolutely do want a good relationship between the hero and the main actress.

Two huge lost opportunities because of no actor personality multicollinearity.

Demo program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols. (My lame blog editor chokes on symbols).

using System;
using System.IO;
using System.Collections.Generic;

namespace VarianceInflationFactor
{
  internal class VarianceInflationFactorProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin variance inflation " +
        "factor (VIF) demo using C# ");

      Console.WriteLine("\nLoading synthetic (20)" +
        " normal data ");
      string file1 = 
        "..\\..\\..\\Data\\synthetic_train_20.txt";
      double[][] train1 = MatLoad(file1, 
        new int[] { 0, 1, 2, 3, 4 }, ',', "#");
      Console.WriteLine("\nFirst three lines: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(train1[i], 4, 9);

      Console.WriteLine("\nBegin VIF analysis ");
      for (int j = 0; j "lt" train1[0].Length; ++j)
      {
        double z = VarInfFactor(train1, j);
        Console.WriteLine("col = " +
          j.ToString().PadLeft(2) + 
          "  |  vif = " + z.ToString("F4"));
      }

      Console.WriteLine("\nLoading synthetic (20)" +
        " highly multicollinear data ");
      Console.WriteLine("(col[2] = 2.0 * col[0] +" +
        " col[1] + noise) ");
      string file2 = 
        "..\\..\\..\\Data\\synthetic_train_20_collinear.txt";
      double[][] train2 = MatLoad(file2,
        new int[] { 0, 1, 2, 3, 4 }, ',', "#");
      Console.WriteLine("\nFirst three lines: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(train2[i], 4, 9);

      Console.WriteLine("\nBegin VIF analysis ");
      for (int j = 0; j "lt" train2[0].Length; ++j)
      {
        double z = VarInfFactor(train2, j);
        Console.WriteLine("col = " +
          j.ToString().PadLeft(2) +
          "  |  vif = " + z.ToString("F4"));
      }
  
      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    // ------------------------------------------------------

    static double VarInfFactor(double[][] data, int col)
    {
      // predict col in data[][] using other columns
      double[][] X = MatWithoutCol(data, col);
      double[] y = MatGetColumn(data, col);
      LinearRegressor model = new LinearRegressor();
      model.Train(X, y);
      //model.TrainLeftPinv(X, y);
      double r2 = model.R2(X, y);  // usually in (0.0, 1.0)
      double vif = 1.0 / (1.0 - r2);  // r2 could be 1
      return vif;
    }

    // ------------------------------------------------------
    // helpers for Main(): MatLoad(), VecShow()
    // ------------------------------------------------------

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      List"lt"double[]"gt" result = 
        new List"lt"double[]"gt"();
      string line = "";
      FileStream ifs = new FileStream(fn, FileMode.Open);
      StreamReader sr = new StreamReader(ifs);
      while ((line = sr.ReadLine()) != null)
      {
        if (line.StartsWith(comment) == true)
          continue;
        string[] tokens = line.Split(sep);
        List"lt"double"gt" lst = new List"lt"double"gt"();
        for (int j = 0; j "lt" usecols.Length; ++j)
          lst.Add(double.Parse(tokens[usecols[j]]));
        double[] row = lst.ToArray();
        result.Add(row);
      }
      sr.Close(); ifs.Close();
      return result.ToArray();
    }

    static void VecShow(double[] vec, int dec, int wid)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      Console.WriteLine("");
    }

    // ------------------------------------------------------
    // helpers for VarInflationFactor()
    // ------------------------------------------------------

    static double[][] MatWithoutCol(double[][] data, int col)
    {
      int nRows = data.Length;
      int nCols = data[0].Length;

      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols-1];

      int k = 0; // into result
      for (int j = 0; j "lt" nCols; ++j)
      {
        if (j == col) continue;
        for (int i = 0; i "lt" nRows; ++i)
          result[i][k] = data[i][j];
        ++k;
      }
      return result;
    }

    static double[] MatGetColumn(double[][] data, int col)
    {
      int nRows = data.Length;
      int nCols = data[0].Length;
      double[] result = new double[nRows];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = data[i][col];
      return result;
    }

    // ------------------------------------------------------

  } // class Program

  // ========================================================

  public class LinearRegressor
  {
    public double[] weights;
    public double bias;
    private Random rnd;

    public LinearRegressor(int seed = 1)
    {
      this.weights = new double[0]; // keep compiler happy
      this.bias = 0;
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------

    public double Predict(double[] x)
    {
      double result = 0.0;
      for (int j = 0; j "lt" x.Length; ++j)
        result += x[j] * this.weights[j];
      result += this.bias;
      return result;
    }

    // ------------------------------------------------------

    public int Train(double[][] trainX, double[] trainY,
      double initRate=0.10, int maxEpochs=10000,
      double noChangeTol=0.00001, int consecutiveNoChange=5)
    {
      // exit when dist(old wts, new wts) lt noChangeTol
      // for consecutiveNoChange times
      int n = trainX.Length;
      int dim = trainX[0].Length;
      this.weights = new double[dim];

      // initialize weights and bias
      // technically not necessary
      double low = -0.01; double hi = 0.01;
      for (int i = 0; i "lt" dim; ++i)
        this.weights[i] = (hi - low) *
          this.rnd.NextDouble() + low;
      this.bias = (hi - low) *
          this.rnd.NextDouble() + low;

      int[] indices = new int[n];  // of train data
      for (int i = 0; i "lt" n; ++i)
        indices[i] = i;

      double[] oldWeightsAndB = new double[dim+1];
      double[] newWeightsAndB = new double[dim+1];
      int countNoChange = 0;
      
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        Shuffle(indices, this.rnd);

        double lrnRate = 
          initRate / Math.Pow((double)(epoch+1), 0.25);
        // Console.WriteLine(lrnRate);

        for (int j = 0; j "lt" dim; ++j)
          oldWeightsAndB[j] = this.weights[j];
        oldWeightsAndB[dim] = this.bias;

        for (int i = 0; i "lt" n; ++i) // each train item
        {
          int idx = indices[i];
          double[] x = trainX[idx];
          double predY = this.Predict(x);
          double actualY = trainY[idx];

          // decay weights before update
          // not needed because regularization not useful
          // for (int j = 0; j "lt" dim; ++j) // each weight
          //   this.weights[j] *= (1 - decay);

          // update weights and bias
          for (int j = 0; j "lt" dim; ++j) // each weight
            this.weights[j] -= lrnRate *
              (predY - actualY) * x[j];
          this.bias -= lrnRate * (predY - actualY) * 1;
        }

        //// display progress 5 times: very noisy
        //if (epoch % (int)(maxEpochs / 5) == 0) // progress
        //{
        //  double r2 = this.R2(trainX, trainY);
        //  string s = "";
        //  s += "epoch = " + epoch.ToString().PadLeft(5);
        //  s += "  R2 = " + r2.ToString("F4").PadLeft(8);
        //  Console.WriteLine(s);
        //}

        // check for early-exit after each epoch
        for (int j = 0; j "lt" dim; ++j)
          newWeightsAndB[j] = this.weights[j];
        newWeightsAndB[dim] = this.bias;
        double d = EuclideanDist(oldWeightsAndB,
          newWeightsAndB);

        if (d "lt" noChangeTol)
        {
          ++countNoChange;
          if (countNoChange == consecutiveNoChange)
            return epoch;
        }
        else
         countNoChange = 0; // reset

      } // epoch
      return maxEpochs;

    } // Train

    // ------------------------------------------------------

    private static double EuclideanDist(double[] v1,
      double[] v2)
    {
      int n = v1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += (v1[i] - v2[i]) * (v1[i] - v2[i]);
      return Math.Sqrt(sum);
    }

    // ------------------------------------------------------

    public double R2(double[][] dataX, double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;

      for (int i = 0; i "lt" n; ++i)
        sum += dataY[i];
      double meanY = sum / n;

      double ssRes = 0.0; // sum squared residuals
      double ssTot = 0.0; // sum squared total
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes +=
          (dataY[i] - predY) * (dataY[i] - predY);
        ssTot +=
          (dataY[i] - meanY) * (dataY[i] - meanY);
      }
      return 1.0 - (ssRes / ssTot);
    }

    // ------------------------------------------------------

    private static void Shuffle(int[] indices, Random rnd)
    {
      int n = indices.Length;
      for (int i = 0; i "lt" n; ++i) // one extra pass
      {
        int ri = rnd.Next(i, n);
        int tmp = indices[i];
        indices[i] = indices[ri];
        indices[ri] = tmp;
      }
    }

  } // class LinearRegressor

} // ns

First, normal, no multicollinearity dataset:

# synthetic_train_20.txt
#
-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
-0.8299, -0.9219, -0.6603,  0.7563, -0.8033,  0.7955
 0.0663,  0.3838, -0.3690,  0.3730,  0.6693,  0.3206
-0.9634,  0.5003,  0.9777,  0.4963, -0.4391,  0.7377
-0.1042,  0.8172, -0.4128, -0.4244, -0.7399,  0.4801
-0.9613,  0.3577, -0.5767, -0.4689, -0.0169,  0.6861
-0.7065,  0.1786,  0.3995, -0.7953, -0.1719,  0.5569
 0.3888, -0.1716, -0.9001,  0.0718,  0.3276,  0.2500
 0.1731,  0.8068, -0.7251, -0.7214,  0.6148,  0.3297
-0.2046, -0.6693,  0.8550, -0.3045,  0.5016,  0.2129
 0.2473,  0.5019, -0.3022, -0.4601,  0.7918,  0.2613
-0.1438,  0.9297,  0.3269,  0.2434, -0.7705,  0.5171
 0.1568, -0.1837, -0.5259,  0.8068,  0.1474,  0.3307
-0.9943,  0.2343, -0.3467,  0.0541,  0.7719,  0.5581
 0.2467, -0.9684,  0.8589,  0.3818,  0.9946,  0.1092
-0.6553, -0.7257,  0.8652,  0.3936, -0.8680,  0.7018
 0.8460,  0.4230, -0.7515, -0.9602, -0.9476,  0.1996

Second, highly multicollinear, dataset:

# synthetic_train_20_collinear.txt
# col [2] = 2*[0] + [1] + rand(0.001)
#
-0.1660,  0.4406,  0.1096, -0.3953, -0.7065, 0.4840
 0.0776, -0.1616, -0.0045, -0.5911,  0.7562, 0.1568
-0.9452,  0.3409, -1.5482,  0.1174, -0.7192, 0.8054
 0.9365, -0.3732,  1.5016,  0.7528,  0.7892, 0.1345
-0.8299, -0.9219, -2.5800,  0.7563, -0.8033, 0.7955
 0.0663,  0.3838,  0.5179,  0.3730,  0.6693, 0.3206
-0.9634,  0.5003, -1.4245,  0.4963, -0.4391, 0.7377
-0.1042,  0.8172,  0.6100, -0.4244, -0.7399, 0.4801
-0.9613,  0.3577, -1.5636, -0.4689, -0.0169, 0.6861
-0.7065,  0.1786, -1.2325, -0.7953, -0.1719, 0.5569
 0.3888, -0.1716,  0.6073,  0.0718,  0.3276, 0.2500
 0.1731,  0.8068,  1.1544, -0.7214,  0.6148, 0.3297
-0.2046, -0.6693, -1.0770, -0.3045,  0.5016, 0.2129
 0.2473,  0.5019,  0.9980, -0.4601,  0.7918, 0.2613
-0.1438,  0.9297,  0.6435,  0.2434, -0.7705, 0.5171
 0.1568, -0.1837,  0.1313,  0.8068,  0.1474, 0.3307
-0.9943,  0.2343, -1.7528,  0.0541,  0.7719, 0.5581
 0.2467, -0.9684, -0.4732,  0.3818,  0.9946, 0.1092
-0.6553, -0.7257, -2.0345,  0.3936, -0.8680, 0.7018
 0.8460,  0.4230,  2.1166, -0.9602, -0.9476, 0.1996