James D. McCaffrey

The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might to predict the maximum safe loan amount to a customer based on customer age, current savings balance, and so on.

There are roughly a dozen major regression techniques, and each technique has several variations. Common techniques are linear regression, linear ridge regression, gradient boosting decision tree regression, and neural network regression. Each technique has pros and cons.

I use the non-standard term bistratum (“two layers”, borrowed from biology) to indicate a neural network with exactly two hidden layers.

I put together a demo using C#. I used one of my standard synthetic datasets. The data 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 first four fields on each line are predictors (often called features). The last value is the target to predict. There are 200 training items and 40 test items.

The output of the demo program is:

Neural network bistratum regression using C#

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 5-(10-10)-1 tanh() identity() neural
 network regressor

Setting lrnRate = 0.0600
Setting maxEpochs = 30000
Setting decay = 0.00000100

Starting training
epoch:      0  MSE = 0.0354  acc = 0.0600
epoch:   6000  MSE = 0.0001  acc = 0.8100
epoch:  12000  MSE = 0.0001  acc = 0.7900
epoch:  18000  MSE = 0.0001  acc = 0.8100
epoch:  24000  MSE = 0.0001  acc = 0.8000
Done

Evaluating model

Accuracy (5%) on train data = 0.8100
Accuracy (5%) on test data  = 0.8000

MSE on train data = 0.0001
MSE on test data = 0.0002

Predicting y for train[0]
Predicted y = 0.4870

End demo

The demo program creates a 5-(10-10)-1 tanh() identity() neural network regression model. This means there are 5 input nodes, two hidden layers with 10 nodes each, and 1 output node. The hidden nodes uses tanh() activation, and the output node uses identity() activation.

The demo program trains the bistratum network using plain vanilla online (aka batch-1) stochastic gradient descent.

The lrnRate (learning rate) influences how quicky the model weights and biases are updated. The maxEpochs controls how many passes through the training data occur. The decay is a weight decay that discourages model weights from becoming extrememly large. When used with SGD, weight decay is mathematically equivalent to L2 regularization.

The model prediction accuracy on the training data is 81% (162 out of 200 correct), and 80% on the test data (32 out of 40 correct). A prediction is scored correct if it is withing 5% of the true target value.

I’m a big fan of 1950s science fiction movies. One of my favorite scenes in many of these old movies is “the monster in the window”.

Right: In “Tarantula” (1955), a biologist in a remote desert lab experiments with ways to grow animals to giant size. What could possibly . . . OK, how about a 50-foot hungry tarantula spider?

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

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

namespace NeuralNetworkBistratumRegression
{
  internal class NeuralNetworkBistratumRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nNeural network bistratum " +
        "regression using C# ");

      // 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 NN
      Console.WriteLine("\nCreating 5-(10-10)-1 tanh()" +
        " identity() neural network regressor ");
      NeuralNetworkBistratumRegressor nn =
         new NeuralNetworkBistratumRegressor(5, 10, 10, 1);

      // 3. train NN
      double lrnRate = 0.06;
      int maxEpochs = 30_000;
      double decay = 0.000001; // .81  .80  .0001 .0002

      Console.WriteLine("\nSetting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpochs = " +
        maxEpochs);
      Console.WriteLine("Setting decay = " +
        decay.ToString("F8"));
      Console.WriteLine("\nStarting training ");
      nn.Train(trainX, trainY, lrnRate, maxEpochs, decay);
      Console.WriteLine("Done ");

      // TODO: programmatically analyze weights and
      // biases for extreme values

      // 4. evaluate trained model
      Console.WriteLine("\nEvaluating model ");
      double trainAcc = nn.Accuracy(trainX, trainY, 0.05);
      double testAcc = nn.Accuracy(testX, testY, 0.05);
      Console.WriteLine("\nAccuracy (5%) on train data = " +
        trainAcc.ToString("F4"));
      Console.WriteLine("Accuracy (5%) on test data  = " +
        testAcc.ToString("F4"));

      double trainMSE = nn.MSE(trainX, trainY);
      double testMSE = nn.MSE(testX, testY);
      Console.WriteLine("\nMSE on train data = " +
        trainMSE.ToString("F4"));
      Console.WriteLine("MSE on test data = " +
        testMSE.ToString("F4"));

      // 5. use model
      Console.WriteLine("\nPredicting y for train[0] ");
      double[] x = trainX[0];
      double predY = nn.Predict(x);
      Console.WriteLine("Predicted y = " +
        predY.ToString("F4"));

      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[][] M)
    {
      int nRows = M.Length;
      int nCols = M[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++] = M[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 NeuralNetworkBistratumRegressor
  {
    // two hidden layers
    public int numInput;
    public int numHiddenA;
    public int numHiddenB;
    public int numOutput;

    public double[] iNodes;  // input nodes
    public double[] aNodes;
    public double[] bNodes;
    public double[] oNodes;  // output nodes

    public double[][] iaWeights; // input to hidden A
    public double[][] abWeights; // hidden A to hidden B
    public double[][] boWeights; // hidden B to output

    public double[] aBiases;
    public double[] bBiases;
    public double[] oBiases;

    private Random rnd; // wt init and train shuffle

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

    public NeuralNetworkBistratumRegressor(int numInput,
      int numHiddenA, int numHiddenB, int numOutput,
      int seed = 0)
    {
      this.numInput = numInput;
      this.numHiddenA = numHiddenA;
      this.numHiddenB = numHiddenB;
      this.numOutput = numOutput; // 1 for regression

      this.iNodes = new double[numInput];
      this.iaWeights = MatMake(numInput, numHiddenA);
      this.aBiases = new double[numHiddenA];
      this.aNodes = new double[numHiddenA];

      this.abWeights = MatMake(numHiddenA, numHiddenB);
      this.bBiases = new double[numHiddenB];
      this.bNodes = new double[numHiddenB];

      this.boWeights = MatMake(numHiddenB, numOutput);
      this.oBiases = new double[numOutput];
      this.oNodes = new double[numOutput];

      this.rnd = new Random(seed);
    }

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

    private 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 double Predict(double[] x)
    {
      // copy input into iNodes
      for (int i = 0; i "lt" this.numInput; ++i)
        this.iNodes[i] = x[i];

      // compute hidden layer A node values
      for (int j = 0; j "lt" this.numHiddenA; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numInput; ++i)
          sum += this.iNodes[i] * this.iaWeights[i][j];
        sum += this.aBiases[j];
        this.aNodes[j] = HyperTan(sum);
      }

      // compute hidden layer B node values
      for (int j = 0; j "lt" this.numHiddenB; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numHiddenA; ++i)
          sum += this.aNodes[i] * this.abWeights[i][j];
        sum += this.bBiases[j];
        this.bNodes[j] = HyperTan(sum);
      }

      // compute output layer node
      for (int j = 0; j "lt" this.numOutput; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numHiddenB; ++i)
          sum += this.bNodes[i] * this.boWeights[i][j];
        sum += this.oBiases[j];
        this.oNodes[j] = Identity(sum);
      }

      return this.oNodes[0];  // a single value
    }

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

    private static double HyperTan(double x)
    {
      if (x "lt" -6.0) return -1.0;
      else if (x "gt" 6.0) return 1.0;
      else return Math.Tanh(x);
    }

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

    private static double Identity(double x)
    {
      return x;
    }

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

    public void Train(double[][] trainX, double[] trainY,
      double lrnRate, int maxEpochs, double decay)
    {
      // 'online' aka batch-1 SGD training (not batch)
      // init weights, leave biases at 0
      double lo = -0.01; double hi = 0.01;

      for (int i = 0; i "lt" this.numInput; ++i)
        for (int j = 0; j "lt" this.numHiddenA; ++j)
          this.iaWeights[i][j] =
            (hi - lo) * this.rnd.NextDouble() + lo;

      for (int i = 0; i "lt" this.numHiddenA; ++i)
        for (int j = 0; j "lt" this.numHiddenB; ++j)
          this.abWeights[i][j] =
            (hi - lo) * this.rnd.NextDouble() + lo;

      for (int i = 0; i "lt" this.numHiddenB; ++i)
        for (int j = 0; j "lt" this.numOutput; ++j)
          this.boWeights[i][j] =
            (hi - lo) * this.rnd.NextDouble() + lo;

      // each weight and bias has a gradient
      double[][] boGrads = MatMake(this.numHiddenB, this.numOutput);
      double[][] abGrads = MatMake(this.numHiddenA, this.numHiddenB);
      double[][] iaGrads = MatMake(this.numInput, this.numHiddenA);

      double[] oBiasGrads = new double[this.numOutput];
      double[] bBiasGrads = new double[this.numHiddenB];
      double[] aBiasGrads = new double[this.numHiddenA];

      // each output and hidden node has a 'signal'
      //  which is gradient without associated input
      //  (lower case delta in Wikipedia)
      double[] oSignals = new double[this.numOutput];
      double[] bSignals = new double[this.numHiddenB];
      double[] aSignals = new double[this.numHiddenA];

      // to shuffle order of items when training
      int[] indices = new int[trainX.Length];
      for (int i = 0; i "lt" indices.Length; ++i)
        indices[i] = i;

      int freq = (int)(maxEpochs / 5); // progress 5 times
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        this.Shuffle(indices); // process in scrambled order
        for (int ii = 0; ii "lt" trainX.Length; ++ii)
        {
          int idx = indices[ii];
          double[] x = trainX[idx];
          double actualY = trainY[idx];
          double predY = this.Predict(x);

          // output node signals depends on target values
          for (int k = 0; k "lt" this.numOutput; ++k)
          {
            double error = predY - actualY;  // standard
            double derivative = 1.0; // identity activation
            oSignals[k] = error * derivative;
          }

          // signal for B nodes depends on output signals
          for (int j = 0; j "lt" this.numHiddenB; ++j)
          {
            // derivative for tanh() activation
            double derivative =
              (1 + this.bNodes[j]) * (1 - this.bNodes[j]);
            double sum = 0.0;
            for (int k = 0; k "lt" this.numOutput; ++k)
              sum += oSignals[k] * this.boWeights[j][k];
            bSignals[j] = derivative * sum;
          }

          // signal for A nodes, depends on B signals
          for (int j = 0; j "lt" this.numHiddenA; ++j)
          {
            // derivative for tanh() activation
            double derivative =
              (1 + this.aNodes[j]) * (1 - this.aNodes[j]);
            double sum = 0.0;
            for (int k = 0; k "lt" this.numHiddenB; ++k)
              sum += bSignals[k] * this.abWeights[j][k];
            aSignals[j] = derivative * sum;
          }

          // at this point, all signals have been computed
          // use signals to calculate gradients left-to-right

          for (int i = 0; i "lt" this.numInput; ++i)
            for (int j = 0; j "lt" this.numHiddenA; ++j)
              iaGrads[i][j] = this.iNodes[i] * aSignals[j];

          for (int i = 0; i "lt" this.numHiddenA; ++i)
            for (int j = 0; j "lt" this.numHiddenB; ++j)
              abGrads[i][j] = this.aNodes[i] * bSignals[j];

          for (int i = 0; i "lt" this.numHiddenB; ++i)
            for (int j = 0; j "lt" this.numOutput; ++j)
              boGrads[i][j] = this.bNodes[i] * oSignals[j];

          // compute bias gradients
          for (int j = 0; j "lt" this.numHiddenA; ++j)
            aBiasGrads[j] = 1.0 * aSignals[j];
          for (int j = 0; j "lt" this.numHiddenB; ++j)
            bBiasGrads[j] = 1.0 * bSignals[j];
          for (int j = 0; j "lt" this.numOutput; ++j)
            oBiasGrads[j] = 1.0 * oSignals[j];

          // before update, decay all weights
          // wt decay is (almost) L2 regularization
          // when used with plain vanilla SGD
          for (int i = 0; i "lt" this.numInput; ++i)
            for (int j = 0; j "lt" this.numHiddenA; ++j)
              this.iaWeights[i][j] *= (1 - decay);

          for (int i = 0; i "lt" this.numHiddenA; ++i)
            for (int j = 0; j "lt" this.numHiddenB; ++j)
              this.abWeights[i][j] *= (1 - decay);

          for (int i = 0; i "lt" this.numHiddenB; ++i)
            for (int j = 0; j "lt" this.numOutput; ++j)
              this.boWeights[i][j] *= (1 - decay);

          // use gradients to update all weights

          for (int i = 0; i "lt" this.numInput; ++i)
            for (int j = 0; j "lt" this.numHiddenA; ++j)
              this.iaWeights[i][j] -= iaGrads[i][j] * lrnRate;

          for (int i = 0; i "lt" this.numHiddenA; ++i)
            for (int j = 0; j "lt" this.numHiddenB; ++j)
              this.abWeights[i][j] -= abGrads[i][j] * lrnRate;

          for (int i = 0; i "lt" this.numHiddenB; ++i)
            for (int j = 0; j "lt" this.numOutput; ++j)
              this.boWeights[i][j] -= boGrads[i][j] * lrnRate;

          // update all biases

          for (int j = 0; j "lt" this.numHiddenA; ++j)
            this.aBiases[j] -= aBiasGrads[j] * lrnRate;

          for (int j = 0; j "lt" this.numHiddenB; ++j)
            this.bBiases[j] -= bBiasGrads[j] * lrnRate;

          for (int j = 0; j "lt" this.numOutput; ++j)
            this.oBiases[j] -= oBiasGrads[j] * lrnRate;

        } // ii each train item

        // display progress
        if (epoch % freq == 0 && epoch "lt" maxEpochs)
        {
          double mse = this.MSE(trainX, trainY);
          double acc = this.Accuracy(trainX, trainY, 0.05);

          string s1 = "epoch: " + 
            epoch.ToString().PadLeft(6);
          string s2 = "  MSE = " + 
            mse.ToString("F4");
          string s3 = "  acc = " + 
            acc.ToString("F4");
          Console.WriteLine(s1 + s2 + s3);
        }

      } // epoch

      return; // no explicit return values
    }

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

    private void Shuffle(int[] indices)
    {
      // Fisher-Yates algorithm
      for (int i = 0; i "lt" indices.Length; ++i)
      {
        int r = this.rnd.Next(i, indices.Length); // note!
        int tmp = indices[r];
        indices[r] = indices[i];
        indices[i] = tmp;
      }
    }

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

    public double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      // in this version, if target y is 0, only pred equal
      // exactly 0 is scored as a correct prediction
      int n = dataX.Length;
      int nCorrect = 0; int nWrong = 0;
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        double actualY = dataY[i];
        if (Math.Abs(predY - actualY) "lte"
          pctClose * Math.Abs(actualY))
          ++nCorrect;
        else
          ++nWrong;
      }
      return (nCorrect * 1.0) / (nCorrect + nWrong);
    }

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

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

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

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

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

    //  double ssRes = 0.0; // sum squares residula
    //  double ssTot = 0.0; // sum squares 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);
    //  }
    //  if (Math.Abs(ssTot) "lt" 1.0e-12) // avoid div by 0
    //    return 0.0;
    //  else
    //    return 1.0 - (ssRes / ssTot);
    //}

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

    //public void SetWeights(double[] weights)
    //{
    //  // used for a Load function, not-yet-implemented
    //  int ptr = 0;
    //  for (int i = 0; i "lt" this.numInput; ++i)
    //    for (int j = 0; j "lt" this.numHiddenA; ++j)
    //      this.iaWeights[i][j] = weights[ptr++];

    //  for (int i = 0; i "lt" numHiddenA; ++i)
    //    this.aBiases[i] = weights[ptr++];

    //  for (int i = 0; i "lt" this.numHiddenA; ++i)
    //    for (int j = 0; j "lt" this.numHiddenB; ++j)
    //      this.abWeights[i][j] = weights[ptr++];

    //  for (int i = 0; i "lt" this.numHiddenB; ++i)
    //    this.bBiases[i] = weights[ptr++];

    //  for (int i = 0; i "lt" this.numHiddenB; ++i)
    //    for (int j = 0; j "lt" this.numOutput; ++j)
    //      this.boWeights[i][j] = weights[ptr++];

    //  for (int i = 0; i "lt" this.numOutput; ++i)
    //    this.oBiases[i] = weights[ptr++];
    //}

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

  } // class NeuralNetworkBistratumRegressor

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

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

Technically, a deep neural network is one that has two or more hidden layers of processing nodes. But deep neural networks for image processing or natural language processing typically have dozens or hundreds of hidden layers.

A special case is a neural network regressor that has exactly two hidden layers. My colleagues and I call such a system a neural network bistratum regressor. Bistratum means two layers in Latin.

In theory, a neural network with a single hidden layer can compute anything a neural network with two ot more layers can, but in practice, two hidden layers often works better.

I put together a demo to illustrate the neural network bistratum regressor input-output mechanism. The output of the demo is:

Neural network bistratum regression IO using C#

Creating 2-3-3-1 tanh() identity() neural network regressor
Done

Setting weights and biases to 0.01 to 0.25
Done

Computing output for x = [1.0, 2.0]
y = 0.4418

Displaying hidden node values:

A: 0.1586 0.1974 0.2355

B: 0.2629 0.2776 0.2923

End demo

The demo network has two input nodes, two hidden layers with three nodes each, and a single output node. The network uses tanh() activation on the hidden layer nodes, and identity() activation on the output node. The network is initialized with weights and biases values from 0.01 to 0.25. The network is fed an input of [1.0, 2.0]. The output is y = 0.4418.

This diagram explains the calculations:

Click to enlarge

The main challenge when using a neural network in finding good values for the weights and biases so that computed outputs match known correct outputs in a set of training data. The most common way to do this is by using one of many variations of stochastic gradient descent.

I’m a big fan of old science fiction movies from the 1950s. Jets had been developed only a few years earlier. Jet bombers made a couple of notable (to me anyway) appearances in two of my favorite movies of the decade.

Bottom Row: In “Beginning of the End” (1957), experiments with radiation to grow giant food result in giant locusts/grasshoppers. Hundreds of thousands of the menaces threaten to destroy Chicago and the move on. As a last resort, the military plans to use a B-36 to drop a nuclear bomb on Chicago. Scientists come up with a clever plan to use a synthetic mating call to lure the insects into Lake Michigan, and humanity is saved.

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

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

namespace NeuralNetworkBistratumRegression
{
  internal class NeuralNetworkBistratumRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nNeural network bistratum " +
        "regression IO using C# ");

      Console.WriteLine("\nCreating 2-3-3-1 tanh()" +
        " identity() neural network regressor ");
      
      NeuralNetworkBistratumRegressor nn =
         new NeuralNetworkBistratumRegressor(2, 3, 3, 1);
      Console.WriteLine("Done ");

      Console.WriteLine("\nSetting weights and biases to" +
        " 0.01 to 0.25 ");
      double[] wts = new double[] {.01, .02, .03, .04, .05,
        .06, .07, .08, .09, .10, .11, .12, .13, .14, .15,
        .16, .17, .18, .19, .20, .21, .22, .23, .24, .25 };
      nn.SetWeights(wts);
      Console.WriteLine("Done ");

      Console.WriteLine("\nComputing output" +
        " for x = [1.0, 2.0] ");
      double y = nn.Predict(new double[] { 1, 2 });
      Console.WriteLine("y = " + y.ToString("F4"));

      Console.WriteLine("\nDisplaying hidden node values: ");
      Console.WriteLine("\nA: " + 
        nn.aNodes[0].ToString("F4") + " " +
        nn.aNodes[1].ToString("F4") + " " +
        nn.aNodes[2].ToString("F4"));
      Console.WriteLine("\nB: " + 
        nn.bNodes[0].ToString("F4") + " " +
        nn.bNodes[1].ToString("F4") + " " +
        nn.bNodes[2].ToString("F4"));

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();

    } // Main

  } // class Program

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

  public class NeuralNetworkBistratumRegressor
  {
    // two hidden layers
    public int numInput;
    public int numHiddenA;
    public int numHiddenB;
    public int numOutput;

    public double[] iNodes;  // input nodes
    public double[] aNodes;
    public double[] bNodes;
    public double[] oNodes;  // output nodes

    public double[][] iaWeights; // input to hidden A
    public double[][] abWeights; // hidden A to hidden B
    public double[][] boWeights; // hidden B to output

    public double[] aBiases;
    public double[] bBiases;
    public double[] oBiases;

    private Random rnd;

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

    public NeuralNetworkBistratumRegressor(int numInput,
      int numHiddenA, int numHiddenB, int numOutput,
      int seed = 0)
    {
      this.numInput = numInput;
      this.numHiddenA = numHiddenA;
      this.numHiddenB = numHiddenB;
      this.numOutput = numOutput; // 1 for regression

      this.iNodes = new double[numInput];
      this.iaWeights = MatMake(numInput, numHiddenA);
      this.aBiases = new double[numHiddenA];
      this.aNodes = new double[numHiddenA];

      this.abWeights = MatMake(numHiddenA, numHiddenB);
      this.bBiases = new double[numHiddenB];
      this.bNodes = new double[numHiddenB];

      this.boWeights = MatMake(numHiddenB, numOutput);
      this.oBiases = new double[numOutput];
      this.oNodes = new double[numOutput];

      this.rnd = new Random(seed);
    }

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

    private 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 double Predict(double[] x)
    {
      // copy input into iNodes
      for (int i = 0; i "lt" this.numInput; ++i)
        this.iNodes[i] = x[i];

      // compute hidden layer A
      for (int j = 0; j "lt" numHiddenA; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numInput; ++i)
          sum += this.iNodes[i] * this.iaWeights[i][j];
        sum += this.aBiases[j];
        this.aNodes[j] = HyperTan(sum);
      }

      // compute hidden layer B
      for (int j = 0; j "lt" this.numHiddenB; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numHiddenA; ++i)
          sum += this.aNodes[i] * this.abWeights[i][j];
        sum += this.bBiases[j];
        this.bNodes[j] = HyperTan(sum);
      }

      // compute output layer node
      for (int j = 0; j "lt" this.numOutput; ++j)
      {
        double sum = 0.0;
        for (int i = 0; i "lt" this.numHiddenB; ++i)
          sum += this.bNodes[i] * this.boWeights[i][j];
        sum += this.oBiases[j];
        this.oNodes[j] = Identity(sum);
      }

      return this.oNodes[0];  // a single value
    }

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

    private static double HyperTan(double x)
    {
      if (x "lt" -6.0) return -1.0;
      else if (x "gt" 6.0) return 1.0;
      else return Math.Tanh(x);
    }

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

    private static double Identity(double x)
    {
      return x;
    }


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

    public void SetWeights(double[] weights)
    {
      int ptr = 0;
      for (int i = 0; i "lt" this.numInput; ++i)
        for (int j = 0; j "lt" this.numHiddenA; ++j)
          this.iaWeights[i][j] = weights[ptr++];

      for (int i = 0; i "lt" numHiddenA; ++i)
        this.aBiases[i] = weights[ptr++];

      for (int i = 0; i "lt" this.numHiddenA; ++i)
        for (int j = 0; j "lt" this.numHiddenB; ++j)
          this.abWeights[i][j] = weights[ptr++];

      for (int i = 0; i "lt" this.numHiddenB; ++i)
        this.bBiases[i] = weights[ptr++];

      for (int i = 0; i "lt" this.numHiddenB; ++i)
        for (int j = 0; j "lt" this.numOutput; ++j)
          this.boWeights[i][j] = weights[ptr++];

      for (int i = 0; i "lt" this.numOutput; ++i)
        this.oBiases[i] = weights[ptr++];
    }

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

  } // class NeuralNetworkBistratumRegressor

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

} // ns

I fed my decision tree regression system, implemented using from-scratch Python with NumPy, 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, using from-scratch Python and NumPy.

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:

Setting max_depth = 3
Setting min_samples = 2
Setting min_leaf = 18
Using default n_split_cols = -1 (all)
Setting save_rows = True

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 associated with node [11]:
[  0   7  14  24  69  87  88 119 121 123 133 136 138 141
 162 186 191 195]

Accuracy train (within 0.10): 0.3750
Accuracy test (within 0.10): 0.4750

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, and the row count associated with the node.

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.

I added an explain() method that shows how a prediction is made:

x = [-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

The explain output has a trailing "AND" with no condition, but I was too lazy to strip it away.

As a sanity check, I ran the synthetic dataset to the scikit-learn DecisionTreeRegressor module and got the same results:

Using scikit with same tree parameters:

Accuracy train (within 0.10): 0.3750
Accuracy test (within 0.10): 0.4750

MSE train: 0.0048
MSE test: 0.0054

x = [-0.1660  0.4406 -0.9998 -0.3953 -0.7065]
Predicted y = 0.4101

Good fun.

I've always been fascinated by models of all kinds. I spent a good portion of my adult career looking at mathematical models. When I was a young man, I was fascinated by model trains. Here's a beautiful HOn30 narrow gauge coal mine.

Demo program. Long, very complex. Replace "lt" (less than), "gt", "lte", "gte" with Boolean operator symbols. (My blog editor chokes on symbols).

# decision_tree_regression_scratch.py

# explicit storage in a list (no pointers/references)
# create iteratively (no stack or recursion)
# performant version at expense of clarity

import numpy as np

class MyDecisionTreeRegressor:  # avoid scikit name

  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 explain(self, x):
    # x is a vector
    curr_idx = 0
    last_valid_value = 0.0
    s = "\nIF \n"
        
    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

      s += "column " + str(curr_node.col_idx) + " "

      if x[curr_node.col_idx] "lte" curr_node.thresh:
        s += " "lte" " + ("%8.4f " % curr_node.thresh)
        s += " AND \n"
        curr_idx = curr_node.left
      else:
        s += " "gt"  " + ("%8.4f " % curr_node.thresh)
        s += " AND \n"
        curr_idx = curr_node.right

    # consider strip away trailing "AND" here . .
    if curr_node is None: nid = -1
    else: nid = curr_node.id
    s += "THEN node [" + str(nid) + "] predicted = "
    s += "%0.4f " % last_valid_value
    print(s)

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

  def display(self):
    for i in range(len(self.tree)):
      n = self.tree[i]
      if n is None: continue

      s1 = "ID %3d " % n.id + " | "
      s2 = "sc %3d " % n.col_idx + " | "
      s3 = "sv %8.4f " % n.thresh + " | "
      s4 = "L %3d " % n.left + " | "
      s5 = "R %3d " % n.right + " | "
      s6 = "py %8.4f " % n.value + " | "
      if n.is_leaf == True: s7 = "leaf T" + " | "
      else: s7 = "leaf F" + " | "
      s8 = "rc %4d " % len(n.rows)
      print(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8)

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

  def best_split(self, rows):
    best_col_idx = -1
    best_thresh = 0.0
    best_var = float('inf')

    n_rows = len(rows)
    if n_rows == 0:
      raise Exception("Empty data in best_split()")

    n_cols = len(self.train_X[0])

    # 1. scramble all column indices (Fisher-Yates)
    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

    # compute safe number of columns to scan, slice out
    if self.n_split_cols == -1:
      n_cols_to_use = n_cols
    else:
      n_cols_to_use = min(self.n_split_cols, n_cols)
    active_cols = col_indices[0:n_cols_to_use]

    # use buffer arrays to avoid memory thrashing in loops
    sorted_rows = np.zeros(n_rows, dtype=np.int32)
    feature_keys = np.zeros(n_rows, dtype=np.float64)

    # calculate sums upfront for entire node population
    total_sum_y = 0.0
    total_sum_sq_y = 0.0
    for i in range(n_rows):
      y_curr = self.train_y[rows[i]]
      total_sum_y += y_curr
      total_sum_sq_y += y_curr * y_curr

    # 2. evaluate each selected column sequentially
    for j in range(len(active_cols)):
      col_idx = active_cols[j]

      # extract the feature values for the active rows
      for i in range(n_rows):
        r = rows[i]
        sorted_rows[i] = r
        feature_keys[i] = self.train_X[r, col_idx]

      # find the index sort order based on feature values
      sort_order = np.argsort(feature_keys)
      
      # apply sort order to buffer arrays
      sorted_features = feature_keys[sort_order]
      sorted_indices = sorted_rows[sort_order]

      left_sum_y = 0.0
      left_sum_sq_y = 0.0

      # 3. loop over the sorted split boundaries
      for i in range(n_rows - 1):
        curr_row_idx = sorted_indices[i]
        y_curr = self.train_y[curr_row_idx]

        left_sum_y += y_curr
        left_sum_sq_y += y_curr * y_curr

        left_count = i + 1
        right_count = n_rows - left_count

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

        current_feature_val = sorted_features[i]
        next_feature_val = sorted_features[i+1]

        # no partitioning identical column values
        if current_feature_val == next_feature_val:
          continue

        right_sum_y = total_sum_y - left_sum_y
        right_sum_sq_y = total_sum_sq_y - left_sum_sq_y

        # fancy math shortcut variance calculation:
        # var = 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 = current_feature_val

    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)

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

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

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

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

    if np.abs(y - y_pred) "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)
  all_preds = model.predict(data_X)
  sum = 0.0
  for i in range(n):
    y = data_y[i]
    y_pred = all_preds[i]
    sum += (y - y_pred) * (y - y_pred)

  return sum / n

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

def main():
  print("\nBegin decision tree regression scratch Python ")

  np.set_printoptions(precision=4, suppress=True,
    floatmode='fixed')
  np.random.seed(0)  # not used this version

  # 1. load data
  print("\nLoading synthetic train (200), test (40) data ")
  train_file = ".\\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
  # . . .

  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: ")
  print(train_X[0:3,:])
  print("\nFirst three y targets: ")
  for i in range(3):
    print("%0.4f" % train_y[i])

  max_depth = 3      # max_depth
  min_samples = 2    # min_samples to consider a split
  min_leaf = 18      # min rows after split
  n_split_cols = -1  # means use all cols
  save_rows = True
  seed = 0

  print("\nSetting max_depth = %1d " % max_depth)
  print("Setting min_samples = %1d " % min_samples)
  print("Setting min_leaf = %1d " % min_leaf)
  print("Using default n_split_cols = -1 (all) ")
  print("Setting save_rows = " + str(save_rows))

  print("\nCreating and training tree ")
  tree = MyDecisionTreeRegressor(max_depth=max_depth, 
    min_samples=min_samples, min_leaf=min_leaf,
    n_split_cols=n_split_cols, save_rows=save_rows,
    seed=seed)
  tree.fit(train_X, train_y)
  print("Done ")

  print("\nTree: ")
  tree.display()

  # assume node [11] exists
  print("\nRows associated with node [11]: ")
  np.set_printoptions(linewidth=60)
  rows = np.array(tree.tree[11].rows, dtype=np.int64)
  print(rows)

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

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

  x = train_X[0]
  print("\nx = ", end=""); print(x)
  y_pred = tree.predict(x.reshape(1,-1))[0]
  print("Predicted y = %0.4f " % y_pred)

  tree.explain(x)

  print("\nEnd demo ")

  print("\n==================== ")

  print("\nUsing scikit with same tree parameters: ")
  from sklearn.tree import DecisionTreeRegressor
  dtr = DecisionTreeRegressor(max_depth=max_depth,
    min_samples_split=min_samples, min_samples_leaf=min_leaf,
    random_state=seed)
  dtr.fit(train_X, train_y)

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

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

  x = train_X[0]
  print("\nx = ", end=""); print(x)
  y_pred = dtr.predict(x.reshape(1,-1))[0]
  print("Predicted y = %0.4f " % y_pred)

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

The goal of a machine learning regression problem is to predict a single numeric value, for example, predicting the bank account balance of a person based on his age, annual income, and so on.

A prediction model can be evaluated in several ways. Four common metrics are: accuracy, mean squared error (MSE), root mean squared error RMSE), and coefficient of determination (aka R2). For accuracy and R2, larger values are better. For MSE and RMSE, smaller values are better.

There are various measures of accuracy, but a typical one is the percentage of correct predictions, where a correct prediction is one that’s within a specified closeness (typically about 5% or 10% or so) to the true target value. Advantage: Easy to interpret. Disadvantage: requires a closeness percentage parameter.

MSE is the average of the squared differences between predicted y and target y values. If the y values have units, such as dollars, MSE has units-squared, such as dollars-squared. RMSE is just the square root of MSE, which, if the y values has units, gives units instead of the awkward units-squared. Advantage: Many regression systems minimize MSE so you get a direct indication of model goodness. Disadvantage: MSE and RMSE values depend on how target values are scaled.

R2 is sort of like accuracy but it doesn’t require a closeness percentage. A better description is that R2 is the proportion of the variance explained by the model. R2 = 1.0 – (SSres / SStot) where SSres = sum(y – y’)^2 and SStot = sum(y – y”)^2. The SSres stands for residual sum of squares. The SStot stands for total sum of squares.

The y is actual target, y’ is predicted target y, and y” is the average of the actual target y values. Yest another way of thinking about R2 is that R2 measures how well the model predicts relative to guessing the average of the target y values. Advantage: Widely used. Disadvantage: Very difficult to interpret meaningfully.

Here’s an implementation of R2 using the C# language (replace “lt” with Boolean less-than symbol). This version is a class method that would be defined inside a class like LinearRegressor or NeuralNetworkRegressor. The implementation assumes there is a class Predict() method.

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

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

  double ssRes = 0.0; // sum squares residula
  double ssTot = 0.0; // sum squares 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);
  }
  if (Math.Abs(ssTot) "lt" 1.0e-12) // avoid div by 0
    return 0.0;
  else
    return 1.0 - (ssRes / ssTot);
}

This method could be called like:

SomeRegressor model = new SomeRegressor();
model.Train(trainX, trainY);
double r2Train = model.R2(trainX, trainY);

An alternative design is to define an external function like:

static double R2(dynamic model, double[][] dataX,
  double[] dataY)
{
  // coefficient of determination
  int n = dataX.Length;
  double sum = 0.0;

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

  double ssRes = 0.0; // sum squares residula
  double ssTot = 0.0; // sum squares total
  for (int i = 0; i "lt" n; ++i) {
    double predY = model.Predict(dataX[i]);
    ssRes += (dataY[i] - predY) * (dataY[i] - predY);
    ssTot += (dataY[i] - meanY) * (dataY[i] - meanY);
  }
  if (Math.Abs(ssTot) "lt" 1.0e-12) // avoid div by 0
    return 0.0;
  else
    return 1.0 - (ssRes / ssTot);
}

The C# “dynamic” keyword allows an object whose type can be determined at runtime. It is analogous to the C# “var” keyword for variables. This external implementation could be called like:

SomeRegressor model = new SomeRegressor();
model.Train(trainX, trainY);
double r2Train = R2(model, trainX, trainY);

One of the main reasons to implement and use an R2 evaluation metric is that almost all of the regression models and classification models in the widely used scikit-learn Python library define a “score” attribute that returns R2. If you implement an R2 metric, you can easily compare non-scikit regression models with scikit regression models.

I’m a big fan of old science fiction movies from the 1950s. Jets had been developed only a few years earlier. Jet bombers made a couple of notable (to me anyway) appearances in two of my favorite movies of the decade.

Bottom Row: In “The Crawling Eye” aka “The Trollenberg Terror” (1958), aliens that look like a cross between a giant eyeball and an octopus, land in the Swiss Alps. The main characters take refuge in a fortified observatory and call in English Electric Canberra bomber to drop napalm on the aliens. The plan succeeds and humanity is saved.

I decided to implement Gradient Boost regression using Blind Trees learners. Bottom line: For my demo dataset, the technique worked about the same as the standard architecture that uses regular decision trees as the learners.

In machine learning, explaining what the problem is, is often more difficult than explaining the solution. So bear with me. But the bottom line is that I tried an experiment that worked, but didn’t provide any significant improvement to the standard approach.

The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might want to predict the maximum safe loan amount to a customer, based on age, account balance, current debt, and so on.

There are many regression techniques. The techniques fall into two main categories: 1.) classical math-based techniques (linear regression, nearest neighbors regression, quadratic regression, kernel ridge regression, neural network regression, and others) and 2.) tree-based techniques (random forest regression, Extra Trees (“extreme random”) regression, AdaBoost regression, Gradient Boost regression, and others).

The tree-based AdaBoost.R2 (adaptive boosting regression, version 2) technique uses a collection of simple decision trees — they’re called the learners or the estimators. Each tree is constructed sequentially, using a different subset of the source training data, with data items that were predicted incorrectly by previous trees are more likely to be included. In this way, each tree gets slightly better. The final prediction is the weighted median of the predictions of the trees.

Although AdaBoost.R2 (often called just AdaBoost) regression almost always uses standard decision trees as the learners, in theory, any kind of simple regression technique can be used. Paradoxically, the base learners need to be weak instead of powerful, but that’s a long and complicated story.

Blind Trees are even weaker than standard decision trees. The tree is completely “blind” to the target y when making its split decisions. It doesn’t calculate target variance, MSE, or target means – it evaluates thresholds purely based on feature space geometry. So Blind Trees trees are much faster than regular decision trees, and they have a kind of built-in regularization that in theory might help AdaBoost regression. (But, alas, didn’t actually help in my experiment).

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 demo are:

Begin AdaBoost.R2 (Blinnd Trees) regression from scratch C#

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

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 nEstimators = 500

Setting lrnRate = 0.5000
Setting tree maxDepth = 8
Setting tree minSamples = 2

Training AdaBoost.R2 model
Done
Created 500 estimators

Accuracy train (within 0.10): 0.9750
Accuracy test (within 0.10): 0.6250

MSE train: 0.0000
MSE test: 0.0018

Predicting for x =
 -0.1660  0.4406 -0.9998 -0.3953 -0.7065
Predicted y = 0.4840

End demo

A very interesting experiment.

AdaBoost regression code is a wrapper around a collection of weak learners. I’m a big fan of science fiction movies of the 1950s. Many of these movies had creatures that were actors wrapped in costumes of some kind. Even though the costumes were weak and not realistic, I still liked many of these movies. Here are two that have an insect theme.

Right: In “The Wasp Woman” (1959), a scientist who works for a cosmetics company develops an anti-aging serum from the royal jelly of wasps. The aging woman owner of the cosmetics firm overdoses on the serum and gets young — but also becomes a murderous wasp-woman. Things don’t end well for her. Like many low-budget sci-fi movies of the 1950s, objectively, this one isn’t very good, but I like it anyway. My grade = C.

Demo program. Very long, quite 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;

namespace AdaBoostWithBlindTrees // AdaBoost.R2 algorithm
{
  internal class AdaBoostWithBlindTreesProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin AdaBoost.R2" +
        " (Blinnd Trees) regression from scratch C# ");

      // 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("\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 model
      int nEstimators = 500;
      double lrnRate = 0.50; // regularizer
      int maxDepth = 8;
      int minSamples = 2;  // .975  .625
      int minLeaf = 1;

      Console.WriteLine("\nSetting nEstimators = " +
        nEstimators);
      Console.WriteLine("\nSetting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting tree maxDepth = " +
        maxDepth);
      Console.WriteLine("Setting tree minSamples = " +
        minSamples);

      Console.WriteLine("\nTraining AdaBoost.R2 model ");
      AdaBoostRegressor model =
        new AdaBoostRegressor(nEstimators, maxDepth,
        minSamples, minLeaf, "linear", lrnRate, seed: 0);
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");
      Console.WriteLine("Created " +
        model.estimators.Count + " estimators ");

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

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

      // 4. use model to make a prediction
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 8);
      double yPred = model.Predict(x);
      Console.WriteLine("Predicted y = " +
        yPred.ToString("F4"));

      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[][] M)
    {
      int nRows = M.Length;
      int nCols = M[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++] = M[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

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

  class AdaBoostRegressor
  {
    public int nEstimators;  // aka nLearners
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public string lossType;
    public double lrnRate;
    public List"lt"BlindTreeRegressor"gt" estimators;
    public List"lt"double"gt" estimatorWeights; // aka alphas
    private Random rnd;

    public AdaBoostRegressor(int nEstimators = 50,
      int maxDepth = 3, int minSamples = 2, int minLeaf=1,
      string lossType = "linear", double lrnRate = 1.0,
      int seed = 0)
    {
      this.nEstimators = nEstimators; // aka learners
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.lossType = lossType;
      this.lrnRate = lrnRate; // not used orig AdaBoost.R2
      this.rnd = new Random(seed);
      this.estimators = new List"lt"BlindTreeRegressor"gt"();
      this.estimatorWeights = new List"lt"double"gt"();
    }

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

    public void Train(double[][] trainX, double[] trainY)
    {
      int nSamples = trainX.Length;

      // 1. initialize uniform sample weights
      double[] weights = new double[nSamples];
      for (int i = 0; i "lt" nSamples; ++i)
        weights[i] = 1.0 / nSamples;

      for (int t = 0; t "lt" this.nEstimators; ++t)
      {
        // normalize weights
        double sumW = 0.0;
        for (int i = 0; i "lt" nSamples; ++i)
          sumW += weights[i];

        double[] wNormed = new double[nSamples];
        for (int i = 0; i "lt" nSamples; ++i)
          wNormed[i] = weights[i] / sumW;

        // draw weighted bootstrap sample using
        // normalized probabilities
        int[] sampleIndices =
          this.MyChoice(nSamples, nSamples, wNormed);

        // get train data subsets
        int nFeatures = trainX[0].Length;
        double[][] subsetX = MatMake(nSamples, nFeatures);
        double[] ySubset = new double[nSamples];

        for (int i = 0; i "lt" nSamples; ++i)
        {
          int idx = sampleIndices[i];
          ySubset[i] = trainY[idx];
          for (int j = 0; j "lt" nFeatures; ++j)
          {
            subsetX[i][j] = trainX[idx][j];
          }
        }

        // 3. train base tree on the bootstrap subset
        int treeSeed = this.rnd.Next(0, 1_000_000);
        BlindTreeRegressor bt =
          new BlindTreeRegressor(
          maxDepth: this.maxDepth,
          minSamples: this.minSamples,
          minLeaf: this.minLeaf,
          numSplitCols: -1,
          saveRows: false,
          seed: treeSeed
        );
        bt.Train(subsetX, ySubset);

        // compute all predictions on full trainX
        double[] preds = new double[nSamples];
        for (int i = 0; i "lt" nSamples; ++i)
          preds[i] = bt.Predict(trainX[i]);

        // absolute errors and max error
        double[] errors = new double[nSamples];
        double maxError = 1e-10; // avoid div by zero
        for (int i = 0; i "lt" nSamples; ++i)
        {
          errors[i] = Math.Abs(preds[i] - trainY[i]);
          if (errors[i] "gt" maxError)
            maxError = errors[i];
        }

        // 5. compute specific loss type
        double[] tLoss = new double[nSamples];
        double eNorm;
        for (int i = 0; i "lt" nSamples; ++i)
        {
          eNorm = errors[i] / maxError; // normalized error
          if (this.lossType == "linear")
            tLoss[i] = eNorm;
          else if (this.lossType == "square")
            tLoss[i] = eNorm * eNorm;
          else
            throw new Exception("unknown loss type ");
        }

        // 6. calculate average weighted error
        double avgError = 0.0;
        for (int i = 0; i "lt" nSamples; ++i)
          avgError += (wNormed[i] * tLoss[i]);

        // if base learner is worse than random guessing,
        // stop boosting
        if (avgError "gte" 0.5)
        {
          if (t == 0)  // first estimator/tree/learner
          {
            this.estimators.Add(bt);
            this.estimatorWeights.Add(1.0e-10);
          }
          break;
        }

        // 7. estimator confidence beta and alpha
        double beta = avgError / (1.0 - avgError);
        if (beta == 0.0) beta = 1.0e-10;

        // moderate estimatorWeights using lrnRate
        double alpha = Math.Log(1.0 / beta);
        this.estimators.Add(bt);
        this.estimatorWeights.Add(this.lrnRate * alpha);

        // 8. update sample weights
        for (int i = 0; i "lt" nSamples; ++i)
        {
          double exp = (1.0 - tLoss[i]);
          double tmp = Math.Pow(beta, exp);
          weights[i] = wNormed[i] * tmp;
        }
      }
    }

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

    private int[] MyChoice(int nItems, int size, double[] p)
    {
      // roulette wheel selection
      // select size ints from [0, nItems) with replacement,
      // using values in vector p as weights

      // default to uniform probability for safety
      //if (p == null)
      //{
      //  p = new double[nItems];
      //  for (int i = 0; i "lt" nItems; ++i)
      //    p[i] = 1.0 / nItems;
      //}

      int[] result = new int[size];

      // compute cumulative distribution function (CDF)
      double[] cdf = new double[nItems];
      double runSum = 0.0;
      for (int i = 0; i "lt" nItems; ++i)
      {
        runSum += p[i];
        cdf[i] = runSum;
      }

      for (int j = 0; j "lt" size; ++j)
      {
        double u = this.rnd.NextDouble();
        int selectedIdx = SearchCdf(cdf, u); // fast binary 
        if (selectedIdx "lt" 0)
          selectedIdx = 0;
        else if (selectedIdx "gte" nItems)
          selectedIdx = nItems - 1;
        result[j] = selectedIdx;
      }

      return result;
    } // MyChoice()

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

    private static int SearchCdf(double[] cdf, double target)
    {
      // binary search to isolate the target interval
      int low = 0;
      int high = cdf.Length - 1;

      while (low "lte" high)
      {
        int mid = low + (high - low) / 2;
        if (cdf[mid] "gte" target)
          high = mid - 1;
        else
          low = mid + 1;
      }

      if (low "gte" cdf.Length) // safety
        return cdf.Length - 1;

      return low;
    }

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

    private 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 double Predict(double[] x)
    {
      int nTrees = this.estimators.Count;
      double[] preds = new double[nTrees];
      double[] modelWts = new double[nTrees];

      for (int t = 0; t "lt" nTrees; ++t)
      {
        preds[t] = this.estimators[t].Predict(x);
        modelWts[t] = this.estimatorWeights[t];
      }

      return WeightedMedian(preds, modelWts);
    }

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

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

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

    public double MSE(double[][] dataX, double[] dataY)
    {
      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;
    }

    // ------------------------------------------------------
    // helper functions for Predict()
    // ------------------------------------------------------

    private static double WeightedMedian(double[] values,
      double[] weights)
    {
      // no interpolation for even n
      // don't assume weights sum to 1.0
      int n = values.Length;
      double sumWts = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sumWts += weights[i];
      double thresh = sumWts / 2;
      int[] sortedIdxs = ArgSort(values);

      double accum = 0.0;
      for (int j = 0; j "lt" n; ++j)
      {
        accum += weights[sortedIdxs[j]];
        if (accum "gte" thresh)
          return values[sortedIdxs[j]];
      }
      return values[sortedIdxs[n - 1]];
    }

    // helper for WeightedMedian()
    private static int[] ArgSort(double[] values)
    {
      int n = values.Length;
      double[] copy = new double[n];
      int[] indices = new int[n];
      for (int i = 0; i "lt" n; ++i)
      {
        copy[i] = values[i];
        indices[i] = i;
      }
      Array.Sort(copy, indices);  // in parallel
      return indices;
    }

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

  } // class AdaBoostRegressor

  // ========================================================
  
  // 1. Target-Blind Selection: The tree is completely "blind"
  // to the target y when making its split decisions.
  // It doesn't calculate target variance, MSE, or target
  // means - it evaluates thresholds purely based on feature
  // space geometry.
  // 2. Feature-Blind Pick: It picks a single feature randomly
  // out of the active set without checking if other features
  // have better predictive power.
  // 3. Guided Blindness: Like a walking stick guiding someone
  // in the dark, using the midpoint prevents the "blind" tree
  // from making silly, extreme splits at the outer edges of
  // data. It moves blindly, but safely through the center of
  // the data space.

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

  public class BlindTreeRegressor
  {
    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 BlindTreeRegressor(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()
    // 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; 
          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;
    }

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

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

      if (nRows == 0)
        throw new Exception("Empty data in BestSplit()");

      // 1. Randomly pick ONE column from active columns
      int nColsToUse = (this.numSplitCols != -1)
          ? Math.Min(this.numSplitCols, nCols)
          : nCols;

      int chosenColIdx = this.rnd.Next(0, nColsToUse);

      // 2. Find min and max for the chosen column
      double minVal = double.MaxValue;
      double maxVal = double.MinValue;
      for (int i = 0; i "lt" nRows; ++i)
      {
        double val = this.trainX[rows[i]][chosenColIdx];
        if (val "lt" minVal) minVal = val;
        if (val "gt" maxVal) maxVal = val;
      }

      if (minVal == maxVal)
        return new double[] { -1.0, 0.0 }; // no split

      // 3. Geometric Midpoint (Fast O(N), no sorting!)
      double baseThresh = (minVal + maxVal) / 2.0;

      // 4. Bias threshold slightly to mean of active rows
      double randThresh = 
        minVal + (maxVal - minVal) * this.rnd.NextDouble();
      double finalThresh = 
        (0.75 * baseThresh) + (0.25 * randThresh);

      // 5. Verify minLeaf condition
      int leftCount = 0;
      int rightCount = 0;
      for (int i = 0; i "lt" nRows; ++i)
      {
        if (this.trainX[rows[i]][chosenColIdx] "lte" 
          finalThresh)
          leftCount++;
        else
          rightCount++;
      }

      if (leftCount "lt" this.minLeaf || 
        rightCount "lt" this.minLeaf)
        return new double[] { -1.0, 0.0 };

      return new double[] { chosenColIdx, finalThresh };
    }

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

    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 BlindTreeRegressor

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

} // 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 wrote an article titled “Quadratic Regression with QR-Householder OLS Solve Training Using C#” in the September 2026 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2026/09/01/quadratic-regression-with-qr-householder-ols-solve-training-using-csharp.aspx.

The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might want to predict the maximum loan amount for a customer based on age, annual salary, bank account balance, and so on. There are approximately a dozen common regression techniques.

Quadratic regression extends basic linear regression. The form of a quadratic regression model is y’ = (w0 * x0) + . . + (wn * xn) + (wj * x0 * x0) + . . + (wk * x0 * x1) + . . . + b. There are derived predictors that are the square of each original predictor, and interaction terms that are the multiplication product of all possible pairs of original predictors.

I implemented a demo using the C# language. The output of the demo is:

Begin C# quadratic regression with direct 
 QR-Householder OLS solver

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 quadratic regression model

Starting direct QR solve training with L2
Setting L2 lamda = 1.0000
Done

Model base weights:
 -0.2590  0.0351 -0.0424  0.0335 -0.1111

Model quadratic weights:
  0.0622  0.0189  0.0054  0.0029  0.0227

Model interaction weights:
  0.0037  0.0245  0.0073  0.1028  0.0007 -0.0097
  0.0347  0.0074 -0.0538  0.0013

Model bias/intercept:   0.3241

Evaluating model
Accuracy train (within 0.10) = 0.8750
Accuracy test (within 0.10) = 0.9500

MSE train = 0.0003
MSE test = 0.0004


Predicting for x =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065

Predicted y = 0.4850

End demo

For my demo, I used one of my standard synthetic datasets. All of the predictor values are between -1 and +1. When using quadratic regression, technically, it’s not necessary to normalize/scale your data. But normalizing is strongly recommended. There are 200 training items and 40 test items.

The demo uses L2 regularization to discourage model overfitting, when the model fits the training data too well and new, previously unseen data is predicted poorly. Technically L2 is not required but is strongly recommended because the process also conditions the data matrix and discourages failure due to arithmetic overflow or underflow.

The trained model predicts the training data with 87.50% accuracy (175 out of 200 correct) and the test data with 95.00% accuracy (38 out of 40 correct). A prediction is scored correct if it’s within 10% of the true target value.

Quadratic regression is most often used with data that has strictly numeric predictor variables. It is possible to use the technique with categorical data, but the details are very tricky and problematic. Many of my colleagues avoid using quadratic regression when the problem scenario has categorical predictor variables.

Quadratic regression is not always effective — if it were, it would be used more often than it is. Compared to basic linear regression, quadratic regression can sometimes provide a big improvement in model prediction accuracy for a relatively small investment in effort, and so it’s usually worth exploring.

I am a big fan of early 1950s science fiction movies. I especially loved scenes where a scientist would explain something using a scale model of some sort.

Right: “Spaceways” (1953) is really a murder mystery and love-triangle story set in the context of a British space program. Here Dr. Smith (facing camera) is a security officer investigating a murder. In the background is a model based on the one designed by rocket scientist Wernher von Braun.