James D. McCaffrey

Bottom line: I compared my from-scratch C# implementation of gradient boosting regression on training data with categorical predictor variables with the scikit-learn GradientBoostingRegression library module and got essentially the same results.

The goal of a machine learning regression problem is to predict a single numeric value. For example, you might want to predict the Income of a person based on their Sex, Age, City, Height, and Psychology (timid, alert, brave). There are many different regression techniques. One of the most common techniques is gradient boosting regression.

Gradient boosting regression (GBR) is a tree-based technique. GBR creates a collection of small decision trees (typically about 100) where each tree learns how to predict training items that were poorly predicted by the previous tree.

I wanted to validate my from-scratch C# implementation of GBR against the Python language scikit-learn library GradientBoostingRegressor module. I generated a set of synthetic training data that looks like:

M, 29, Austin, 68.5, timid, 54600.00
F, 51, Juneau, 60.0, brave, 25700.00
M, 32 ,Boston, 73.5, alert, 74900.00
. . .

The goal is to predict Income (in the last column) from Sex, Age, City, Height (in inches), and Psychology (timid, alert, brave). There are 200 training items and 40 test items.

The Sex variable is encoded as M = 0, F = 1. There are 10 possible City values: Austin, Boston, Casper, Denver, Eugene, Fresno, Gallup, Helena, Irvine, Juneau.

I used sorted-target-mean ordinal encoding, where the city that has the smallest average Income (Boston = $30,600.00) is encoded as 0, the next larger average Income city (Casper = $33,100.00) is encoded as 1, and so on to the largest average income city (Helena = $64,194.74) is 9.

I used one-hot encoding for Psychology: timid = 1,0,0 alert = 0,1,0 brave = 0,0,1.

When using tree-based regression, it’s not necessary to normalize numeric predictor variables, but it doesn’t hurt and can be helpful in some scenarios. I divided all Age values by 100, and all Height values by 100. I normalized the target Income variable by dividing all values by 100,000.

The resulting normalized and encoded data looks like:

0, 0.29, 2, 0.685, 1,0,0, 0.54600
1, 0.51, 8, 0.600, 0,0,1, 0.25700
0, 0.32, 0, 0.735, 0,1,0, 0.74900
. . . 

The key parts of the output of my C# implementation of gradient boosting regression are:

Begin C# Gradient Boosting regression demo
Predict Income from Sex, Age, City Height, Psychology

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

First three raw train X:
M,29,Austin,68.5,timid,54600.00
F,51,Juneau,60.0,brave,25700.00
M,32,Boston,73.5,alert,74900.00

First three train X:
  0.00000  0.29000  2.00000  0.68500  1.00000  0.00000  0.00000
  1.00000  0.51000  8.00000  0.60000  0.00000  0.00000  1.00000
  0.00000  0.32000  0.00000  0.73500  0.00000  1.00000  0.00000

First three train y:
  0.54600
  0.25700
  0.74900

Setting numTrees = 200
Setting maxDepth = 4
Setting minSamples = 2
Setting minLeaf = 1
Setting lrnRate = 0.0500

Creating and training GradientBoostRegression model
Done

Evaluating model

Accuracy train (within 0.10) = 1.0000
Accuracy test (within 0.10) = 0.7750

MSE train = 0.0000
MSE test = 0.0022

R2 train = 0.9993
R2 test = 0.9344

Predicting for x =
   0.0000   0.2900   2.0000   0.6850   1.0000   0.0000   0.0000
Predicted y = 0.5391

Explaining for x =
   0.0000   0.2900   2.0000   0.6850   1.0000   0.0000   0.0000

Initial prediction: 0.4685
t =   0  pred_res =  -0.0342  delta =  -0.0017  pred =  0.4668
t =   1  pred_res =   0.0079  delta =   0.0004  pred =  0.4672
t =   2  pred_res =  -0.0315  delta =  -0.0016  pred =  0.4656
t =   3  pred_res =  -0.0361  delta =  -0.0018  pred =  0.4638
t =   4  pred_res =  -0.0136  delta =  -0.0007  pred =  0.4632
. . .
t = 198  pred_res =  -0.0000  delta =  -0.0000  pred =  0.5391
t = 199  pred_res =   0.0000  delta =   0.0000  pred =  0.5391
Predicted y = 0.5391

End demo

I ran the data through the scikit GradientBoostingRegression module and got essentially the same results:

Begin scikit gradient boost regression demo
Predict Income from Sex, Age, City, Height, Psych

Setting n_estimators = 200
Setting max_depth = 4
Setting min_samples_split = 2
Setting min_samples_leaf = 1
Setting lrn_rate = 0.0500

Training gradient boost regression model
Done

Evaluating model

Acc train (within 0.10): 1.0000
Acc test (within 0.10): 0.8000

MSE train: 0.0000
MSE test: 0.0022

R2 train: 0.9993
R2 test: 0.9347

Predicting for x =
[0.0000 0.2900 2.0000 0.6850 1.0000 0.0000 0.0000]
Predicted y = 0.5391

End demo

A big advantage of my from-scratch C# implementation of GBR is that it was easy for me to add statements that show exactly how a prediction is made.

Machine learning regression models are often beautiful to my eye. But I’m not a good judge of other kinds of models. I’ve noticed that fashion models often look like they walk with an exaggerated limp. Here are three images that were returned by an Internet search for “ugly fashion model on runway”.

Right: OK, there isn’t much to disagree about “ugly fashion model on runway” here. I don’t want to sound mean, but she doesn’t belong on a fashion runway. An airport runway maybe, but not a fashion runway.

Demo C# program. Very long and complex. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my lame blog editor consistently chokes on symbols).

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

namespace GradientBoostRegression
{
  internal class GradientBoostRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# Gradient Boosting" +
        " regression demo ");

      Console.WriteLine("Predict Income from Sex, Age," +
        " City Height, Psychology ");

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

      string testFile =
        "..\\..\\..\\Data\\People\\people_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        [7], ',', "#"));
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three raw train X: ");
      Console.WriteLine("M,29,Austin,68.5,timid,54600.00 ");
      Console.WriteLine("F,51,Juneau,60.0,brave,25700.00 ");
      Console.WriteLine("M,32,Boston,73.5,alert,74900.00 ");

      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 numTrees = 200;
      int maxDepth = 4;
      int minSamples = 2;
      int minLeaf = 1;
      int numSplitCols = -1; // use all columns
      double lrnRate = 0.10;


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

      Console.WriteLine("\nCreating and training" +
        " GradientBoostRegression model ");
      GradientBoostRegressor gbr =
        new GradientBoostRegressor(numTrees, maxDepth,
        minSamples, minLeaf, numSplitCols, lrnRate);
      gbr.Train(trainX, trainY);
      Console.WriteLine("Done ");

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

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

      double r2Train = gbr.R2(trainX, trainY);
      Console.WriteLine("\nR2 train = " +
        r2Train.ToString("F4"));
      double r2Test = gbr.R2(testX, testY);
      Console.WriteLine("R2 test = " +
        r2Test.ToString("F4"));

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

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

    // ------------------------------------------------------
    // helpers for Main()
    // ------------------------------------------------------

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

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

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

  } // class Program

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

  public class GradientBoostRegressor
  {
    public double lrnRate;
    public int nTrees;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"DecisionTreeRegressor"gt" trees;
    public double pred0;  // initial prediction
    private Random rnd;

    public GradientBoostRegressor(int nTrees, int maxDepth,
      int minSamples, int minLeaf, int numSplitCols,
      double lrnRate, int seed=0)
    {
      this.nTrees = nTrees;
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.trees = new List"lt"DecisionTreeRegressor"gt"();
      this.lrnRate = lrnRate;
      this.rnd = new Random(seed);
    }

    public void Train(double[][] trainX, double[] trainY)
    {
      int n = trainX.Length;
      this.pred0 = Mean(trainY);

      double[] preds = new double[n]; //each data item
      for (int i = 0; i "lt" n; ++i)
        preds[i] = this.pred0;

      double[] residuals = new double[n]; // for curr tree
      for (int t = 0; t "lt" this.nTrees; ++t) // each tree
      {
        // double[] residuals = new double[n]; // curr tree
        for (int i = 0; i "lt" n; ++i)
          residuals[i] = trainY[i] - preds[i];

        DecisionTreeRegressor dtr =
          new DecisionTreeRegressor(this.maxDepth,
          this.minSamples, this.minLeaf, this.numSplitCols,
          false, this.rnd.Next(0,1_000_000));
        dtr.Train(trainX, residuals); // predict residuals

        for (int i = 0; i "lt" n; ++i)
        {
          double predResidual = dtr.Predict(trainX[i]);
          preds[i] += this.lrnRate * predResidual;
        }
        this.trees.Add(dtr);
      }
    } // Train

    public double Predict(double[] x, bool verbose = false)
    {
      double result = this.pred0;
      if (verbose == true)
      {
        Console.WriteLine("\nInitial prediction: " +
        result.ToString("F4"));
      }
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;

        if (verbose == true)
        {
          if (t "gte" 0 && t "lte" 4 ||
            t "gte" this.nTrees - 2 &&
            t "lte" this.nTrees - 1)
          {
            Console.Write("t = " + t.ToString().PadLeft(3) +
              "  pred_res = " + predResidual.ToString("F4").
              PadLeft(10));
            Console.Write("  delta = " + delta.ToString("F4").
              PadLeft(10));
            Console.WriteLine("  pred = " +
              result.ToString("F4").PadLeft(7));
          }
          if (t == 5)
            Console.WriteLine(". . . ");
        }
      }
      return result;
    }

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

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

    public 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)
        sum += dataY[i];
      double meanActual = sum / n;

      double ssRes = 0.0;
      double ssTot = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes += (dataY[i] - predY) * (dataY[i] - predY);
        ssTot += (dataY[i] - meanActual) *
          (dataY[i] - meanActual);
      }
      if (Math.Abs(ssTot) "lt" 1.0e-12)
        return 0.0;
      else
        return 1.0 - (ssRes / ssTot);
    }

    private static double Mean(double[] data)
    {
      int n = data.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += data[i];
      return sum / n;
    }

  } // class GradientBoostRegressor

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

  public class DecisionTreeRegressor
  {
    // ChatGPT: "I don't see any fundamental tree-building or
    // indexing bugs. If the inputs are guaranteed to be
    // finite doubles and the parameters are validated 
    // elsewhere, the core regression algorithm is correct."

    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 = -1;
      public double thresh = 0.0;
      public int left = -1;
      public int right = -1;
      public double value = 0.0;
      public bool isLeaf = false;
      public List"lt"int"gt" rows;
    }

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

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

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

    // ------------------------------------------------------
    // ------------------------------------------------------
    // public: ctor(), Train(), Predict()
    // private: BestSplit(), TreeTargetMean()
    // ------------------------------------------------------

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

      int maxParentID = (int)Math.Pow(2, this.maxDepth) - 2;

      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);

      // construct root node
      Node root = new Node()
      {
        id = 0,
        value = this.TreeTargetMean(allRows),
        rows = allRows
      };
      this.tree[0] = root;

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

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

        // 1. check before expensive split calcs
        if (currNode.id "gt" maxParentID ||
            currNode.rows.Count "lt" this.minSamples ||
            rightID "gte" this.tree.Count)
        {
          currNode.isLeaf = true;
          continue;
        }

        // 2. find best split
        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1];

        if (colIdx == -1) // unable to split
        {
          currNode.isLeaf = true;
          continue;
        }

        // 3. apply split properties
        currNode.colIdx = colIdx;  // aka split col
        currNode.thresh = splitVal;

        List"lt"int"gt" leftIdxs = new List"lt"int"gt"();
        List"lt"int"gt" rightIdxs = new List"lt"int"gt"();

        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);
        }

        // 4. create and attach children
        currNode.left = leftID;
        this.tree[leftID] = new Node()
        {
          id = leftID,
          rows = leftIdxs,
          value = this.TreeTargetMean(leftIdxs)
        };

        currNode.right = rightID;
        this.tree[rightID] = new Node()
        {
          id = rightID,
          rows = rightIdxs,
          value = this.TreeTargetMean(rightIdxs)
        };
      } // main loop

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

    } // Train()

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

    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 || 
          currNode.colIdx == -1) 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 bestColIdx = -1;
      double bestThresh = 0.0;
      double bestGain = -1.0;
      int nRows = rows.Count;
      int nCols = this.trainX[0].Length;

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

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

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

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

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

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

        Array.Sort(featureKeys, sortedRows);
        double leftSum = 0.0;

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

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

          if (featureKeys[i] == featureKeys[i + 1])
            continue; // duplicate values
          if (leftCount "lt" this.minLeaf ||
            rightCount "lt" this.minLeaf)
            continue;

          double rightSum = totalSum - leftSum;

          // maximizing Gain = (G_L^2 / N_L) + (G_R^2 / N_R)
          // equivalent to maximizing explicit weighted
          // variance reduction
          double gain = (leftSum * leftSum / leftCount) +
            (rightSum * rightSum / rightCount);

          if (gain "gt" bestGain)
          {
            bestGain = gain;
            bestColIdx = colIdx;
            // bestThresh = 
            //  (featureKeys[i] + featureKeys[i + 1]) / 2.0;
            bestThresh = featureKeys[i] +
              (featureKeys[i + 1] - featureKeys[i]) / 2.0;
          }
        }
      }

      // using a tuple would create version dependency
      double[] result = new double[2];
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;
    }

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

  } // DecisionTreeRegressor

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

  

Demo scikit program.

# scikit_gradient_boost_regression.py
# on custom People dataset

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor

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

def accuracy(model, data_X, data_y, pct_close):
  # assumes model has a predict(X)
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)  # make it a matrix
    y = data_y[i]
    y_pred = model.predict(x)[0]  # predict() expects 2D

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

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

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

  return sum / n

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

def main():
  print("\nBegin scikit gradient boost regression demo ")
  print("Predict Income from Sex, Age, City, Height, Psych ")

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

  # 1. load data
  print("\nLoading People train (200), test (400) data ")
 
  train_file = ".\\Data\\People\\people_train_200.txt"
  test_file = ".\\Data\\People\\people_test_40.txt"
  # sex, age, city, height, psych(3)
  x_cols=[0,1,2,3, 4,5,6]
  target_col = 7  # 1.0000 .6413

  train_X = np.loadtxt(train_file, comments="#",
    usecols=x_cols, delimiter=",",  dtype=np.float64)
  train_y = np.loadtxt(train_file, comments="#",
    usecols=target_col, delimiter=",",  dtype=np.float64)

  test_X = np.loadtxt(test_file, comments="#",
    usecols=x_cols, delimiter=",",  dtype=np.float64)
  test_y = np.loadtxt(test_file, comments="#",
    usecols=target_col, delimiter=",",  dtype=np.float64)
  print("Done ")

  print("\nFirst three raw train X: ")
  print("M,29,Austin,68.5,timid,54600.00 ")
  print("F,51,Juneau,60.0,brave,25700.00 ")
  print("M,32,Boston,73.5,alert,74900.00 ")

  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])

  n_estimators = 200
  max_depth = 4
  min_samples_split = 2
  min_samples_leaf = 1
  lrn_rate = 0.05
  seed = 0

  print("\nSetting n_estimators = " + \
    str(n_estimators))
  print("Setting max_depth = " + \
    str(max_depth))
  print("Setting min_samples_split = " + \
    str(min_samples_split))
  print("Setting min_samples_leaf = " + \
    str(min_samples_leaf))
  print("Setting lrn_rate = %0.4f " % lrn_rate)

  # GradientBoostingRegressor(*, loss='squared_error',
  # learning_rate=0.1, n_estimators=100, subsample=1.0,
  # criterion='friedman_mse', min_samples_split=2,
  # min_samples_leaf=1, min_weight_fraction_leaf=0.0,
  # max_depth=3, min_impurity_decrease=0.0, init=None,
  # random_state=None, max_features=None, alpha=0.9,
  # verbose=0, max_leaf_nodes=None, warm_start=False,
  # validation_fraction=0.1, n_iter_no_change=None,
  # tol=0.0001, ccp_alpha=0.0)

  model = \
    GradientBoostingRegressor(max_depth=max_depth, 
    min_samples_split=min_samples_split,
    n_estimators=n_estimators,
    min_samples_leaf=min_samples_leaf,
    learning_rate=lrn_rate, random_state=seed)

  print("\nTraining gradient boost regression model ")
  model.fit(train_X, train_y)
  print("Done " )

  print("\nEvaluating model ")
  acc_train = accuracy(model, train_X, train_y, 0.10)
  print("\nAcc train (within 0.10): %0.4f " % acc_train)
  acc_test = accuracy(model, test_X, test_y, 0.10)
  print("Acc test (within 0.10): %0.4f " % acc_test)

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

  r2_train = model.score(train_X, train_y)
  print("\nR2 train: %0.4f " % r2_train)
  r2_test = model.score(test_X, test_y)
  print("R2 test: %0.4f " % r2_test)

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

  print("\nEnd demo ")

if __name__ == "__main__":
  main()

Training data:

# people_train_200.txt
#
# M = 0,  F = 1
# city_target_encoding = {
#     "Boston": 30600.00,  0
#     "Casper": 33100.00,  1
#     "Austin": 33961.90,  2
#     "Denver": 35535.00,  3
#     "Eugene": 44090.00,  4
#     "Irvine": 51840.00,  5
#     "Fresno": 52950.00,  6
#     "Gallup": 52960.00,  7
#     "Juneau": 58842.86,  8
#     "Helena": 64194.74,  9
# }
# timid = 1,0,0  alert = 0,1,0  brave = 0,0,1
#
# Sex Age City Height Psych Income
#
0, 0.29, 2, 0.685, 1,0,0, 0.54600
1, 0.51, 8, 0.600, 0,0,1, 0.25700
0, 0.32, 0, 0.735, 0,1,0, 0.74900
0, 0.19, 5, 0.610, 0,0,1, 0.30300
1, 0.45, 4, 0.645, 0,1,0, 0.35600
0, 0.58, 6, 0.720, 1,0,0, 0.52100
1, 0.22, 2, 0.630, 0,1,0, 0.27800
0, 0.64, 8, 0.775, 0,0,1, 0.88200
1, 0.38, 3, 0.660, 1,0,0, 0.32400
0, 0.41, 7, 0.700, 0,1,0, 0.51200
1, 0.27, 1, 0.625, 0,0,1, 0.28900
0, 0.53, 9, 0.740, 1,0,0, 0.62300
1, 0.34, 0, 0.650, 0,1,0, 0.33100
0, 0.60, 5, 0.765, 0,0,1, 0.79800
1, 0.21, 4, 0.605, 1,0,0, 0.22400
0, 0.48, 6, 0.690, 0,1,0, 0.48700
1, 0.62, 2, 0.670, 0,0,1, 0.53600
0, 0.25, 8, 0.715, 1,0,0, 0.45200
1, 0.39, 3, 0.640, 0,1,0, 0.31800
0, 0.55, 7, 0.750, 0,0,1, 0.72100
1, 0.31, 1, 0.615, 1,0,0, 0.24500
0, 0.43, 9, 0.705, 0,1,0, 0.52400
1, 0.50, 0, 0.635, 0,0,1, 0.38900
0, 0.28, 5, 0.730, 0,1,0, 0.58100
1, 0.56, 4, 0.665, 1,0,0, 0.47800
0, 0.20, 6, 0.620, 0,0,1, 0.29100
1, 0.42, 2, 0.655, 0,1,0, 0.37200
0, 0.66, 8, 0.780, 1,0,0, 0.85400
1, 0.36, 3, 0.620, 0,0,1, 0.29800
0, 0.49, 7, 0.725, 0,1,0, 0.59300
1, 0.24, 1, 0.600, 0,1,0, 0.21200
0, 0.57, 9, 0.755, 0,0,1, 0.76400
1, 0.46, 0, 0.640, 1,0,0, 0.36100
0, 0.33, 5, 0.680, 0,1,0, 0.44300
1, 0.61, 4, 0.675, 0,0,1, 0.57900
0, 0.23, 6, 0.710, 1,0,0, 0.42800
1, 0.52, 2, 0.630, 0,1,0, 0.39400
0, 0.65, 8, 0.760, 0,0,1, 0.81200
1, 0.30, 3, 0.610, 1,0,0, 0.23100
0, 0.44, 7, 0.695, 0,1,0, 0.49600
1, 0.59, 1, 0.680, 0,0,1, 0.58300
0, 0.35, 9, 0.720, 1,0,0, 0.49800
1, 0.18, 0, 0.605, 0,1,0, 0.21900
0, 0.51, 5, 0.745, 0,0,1, 0.66700
1, 0.40, 4, 0.650, 1,0,0, 0.33500
0, 0.67, 6, 0.785, 0,1,0, 0.94200
1, 0.26, 2, 0.625, 0,0,1, 0.27400
0, 0.37, 8, 0.690, 0,1,0, 0.45100
1, 0.54, 3, 0.660, 1,0,0, 0.43600
0, 0.30, 7, 0.730, 0,0,1, 0.59800
1, 0.47, 1, 0.645, 0,1,0, 0.38200
0, 0.63, 9, 0.770, 1,0,0, 0.83100
1, 0.28, 0, 0.615, 0,0,1, 0.26300
0, 0.42, 5, 0.700, 0,1,0, 0.48900
1, 0.35, 4, 0.635, 1,0,0, 0.28700
0, 0.59, 6, 0.750, 0,0,1, 0.74600
1, 0.20, 2, 0.610, 0,1,0, 0.22800
0, 0.46, 8, 0.720, 1,0,0, 0.53200
1, 0.58, 3, 0.670, 0,1,0, 0.51400
0, 0.31, 7, 0.715, 0,0,1, 0.52900
1, 0.33, 1, 0.620, 1,0,0, 0.25800
0, 0.54, 9, 0.735, 0,1,0, 0.64500
1, 0.41, 0, 0.655, 0,0,1, 0.36700
0, 0.22, 5, 0.630, 1,0,0, 0.26400
1, 0.67, 4, 0.685, 0,1,0, 0.69100
0, 0.38, 6, 0.695, 0,0,1, 0.47900
1, 0.25, 2, 0.600, 1,0,0, 0.19800
0, 0.61, 8, 0.760, 0,1,0, 0.77400
1, 0.43, 3, 0.640, 0,0,1, 0.35200
0, 0.26, 7, 0.700, 1,0,0, 0.43100
1, 0.48, 1, 0.650, 0,1,0, 0.39700
0, 0.50, 9, 0.725, 0,0,1, 0.60200
1, 0.29, 0, 0.620, 1,0,0, 0.24900
0, 0.68, 5, 0.785, 0,1,0, 0.96800
1, 0.37, 4, 0.630, 0,0,1, 0.30600
0, 0.18, 6, 0.615, 1,0,0, 0.25300
1, 0.65, 2, 0.680, 0,1,0, 0.63200
0, 0.34, 8, 0.685, 0,0,1, 0.43800
1, 0.53, 3, 0.655, 1,0,0, 0.41900
0, 0.45, 7, 0.710, 0,1,0, 0.53700
1, 0.23, 1, 0.610, 0,0,1, 0.23600
0, 0.52, 9, 0.740, 1,0,0, 0.61100
1, 0.32, 0, 0.630, 0,1,0, 0.28400
0, 0.62, 5, 0.770, 0,0,1, 0.82500
1, 0.19, 4, 0.600, 1,0,0, 0.20100
0, 0.40, 6, 0.700, 0,1,0, 0.47300
1, 0.60, 2, 0.675, 0,0,1, 0.56800
0, 0.27, 8, 0.695, 1,0,0, 0.41200
1, 0.44, 3, 0.645, 0,1,0, 0.36900
0, 0.56, 7, 0.755, 0,0,1, 0.75300
1, 0.36, 1, 0.635, 1,0,0, 0.29100
0, 0.47, 9, 0.715, 0,1,0, 0.55800
1, 0.55, 0, 0.660, 0,0,1, 0.46200
0, 0.21, 5, 0.625, 0,1,0, 0.28900
1, 0.63, 4, 0.680, 1,0,0, 0.58400
0, 0.39, 6, 0.690, 0,0,1, 0.46100
1, 0.22, 2, 0.615, 1,0,0, 0.21700
0, 0.64, 8, 0.770, 0,1,0, 0.85900
1, 0.38, 3, 0.635, 0,0,1, 0.33800
0, 0.29, 7, 0.705, 1,0,0, 0.45800
1, 0.51, 1, 0.655, 0,1,0, 0.43100
0, 0.58, 9, 0.760, 0,0,1, 0.78200
1, 0.27, 0, 0.610, 0,1,0, 0.24300
0, 0.36, 5, 0.685, 1,0,0, 0.42100
1, 0.66, 4, 0.685, 0,0,1, 0.67500
0, 0.24, 6, 0.640, 0,1,0, 0.31200
1, 0.49, 2, 0.650, 1,0,0, 0.38400
0, 0.60, 8, 0.755, 0,0,1, 0.76900
1, 0.31, 3, 0.625, 0,1,0, 0.26800
0, 0.43, 7, 0.700, 1,0,0, 0.48200
1, 0.57, 1, 0.665, 0,0,1, 0.51900
0, 0.32, 9, 0.690, 0,1,0, 0.43900
1, 0.20, 0, 0.600, 1,0,0, 0.19200
0, 0.53, 5, 0.750, 0,1,0, 0.68400
1, 0.42, 4, 0.640, 0,0,1, 0.36100
0, 0.65, 6, 0.775, 1,0,0, 0.84700
1, 0.28, 2, 0.620, 0,1,0, 0.26200
0, 0.35, 8, 0.680, 1,0,0, 0.39600
1, 0.50, 3, 0.650, 0,0,1, 0.41800
0, 0.25, 7, 0.690, 0,1,0, 0.41300
1, 0.45, 1, 0.640, 1,0,0, 0.34800
0, 0.59, 9, 0.765, 0,1,0, 0.80100
1, 0.33, 0, 0.625, 0,0,1, 0.30200
0, 0.23, 5, 0.635, 0,1,0, 0.29700
1, 0.64, 4, 0.680, 0,0,1, 0.62100
0, 0.41, 6, 0.705, 1,0,0, 0.47100
1, 0.21, 2, 0.605, 0,0,1, 0.21800
0, 0.67, 8, 0.785, 0,1,0, 0.95600
1, 0.37, 3, 0.630, 1,0,0, 0.29200
0, 0.48, 7, 0.720, 0,0,1, 0.59100
1, 0.26, 1, 0.615, 0,1,0, 0.26100
0, 0.55, 9, 0.745, 1,0,0, 0.65800
1, 0.40, 0, 0.645, 0,1,0, 0.34900
0, 0.19, 5, 0.615, 0,0,1, 0.28400
1, 0.68, 4, 0.690, 1,0,0, 0.70200
0, 0.36, 6, 0.680, 0,1,0, 0.42600
1, 0.29, 2, 0.625, 0,0,1, 0.29300
0, 0.62, 8, 0.765, 1,0,0, 0.79300
1, 0.44, 3, 0.640, 0,1,0, 0.35700
0, 0.28, 7, 0.700, 0,0,1, 0.46900
1, 0.52, 1, 0.660, 1,0,0, 0.42100
0, 0.46, 9, 0.710, 0,1,0, 0.53900
1, 0.24, 0, 0.605, 0,0,1, 0.20700
0, 0.57, 5, 0.755, 1,0,0, 0.75100
1, 0.35, 4, 0.630, 0,1,0, 0.29600
0, 0.66, 6, 0.780, 0,0,1, 0.88900
1, 0.30, 2, 0.620, 0,0,1, 0.29700
0, 0.38, 8, 0.690, 1,0,0, 0.42400
1, 0.56, 3, 0.665, 0,1,0, 0.48300
0, 0.27, 7, 0.695, 0,0,1, 0.44800
1, 0.43, 1, 0.635, 1,0,0, 0.32800
0, 0.61, 9, 0.770, 0,1,0, 0.83600
1, 0.31, 0, 0.620, 0,0,1, 0.28100
0, 0.20, 5, 0.620, 1,0,0, 0.24600
1, 0.65, 4, 0.685, 0,1,0, 0.66400
0, 0.42, 6, 0.700, 0,0,1, 0.51200
1, 0.18, 2, 0.600, 0,1,0, 0.18900
0, 0.63, 8, 0.775, 1,0,0, 0.82900
1, 0.39, 3, 0.635, 0,1,0, 0.32600
0, 0.50, 7, 0.725, 0,0,1, 0.63400
1, 0.25, 1, 0.610, 1,0,0, 0.22300
0, 0.54, 9, 0.740, 0,1,0, 0.65200
1, 0.36, 0, 0.635, 1,0,0, 0.29800
0, 0.31, 5, 0.675, 0,0,1, 0.42900
1, 0.58, 4, 0.670, 0,1,0, 0.52800
0, 0.25, 6, 0.645, 1,0,0, 0.31900
1, 0.47, 2, 0.645, 0,0,1, 0.41100
0, 0.59, 8, 0.750, 0,1,0, 0.73800
1, 0.28, 3, 0.615, 1,0,0, 0.23700
0, 0.40, 7, 0.690, 0,1,0, 0.45600
1, 0.54, 1, 0.650, 0,0,1, 0.46300
0, 0.34, 9, 0.695, 1,0,0, 0.43200
1, 0.23, 0, 0.600, 0,1,0, 0.20800
0, 0.52, 5, 0.740, 0,0,1, 0.66100
1, 0.41, 4, 0.640, 1,0,0, 0.32900
0, 0.68, 6, 0.785, 0,1,0, 0.97800
1, 0.26, 2, 0.610, 0,0,1, 0.23800
0, 0.37, 8, 0.685, 0,1,0, 0.40900
1, 0.50, 3, 0.655, 1,0,0, 0.42700
0, 0.29, 7, 0.700, 0,0,1, 0.46300
1, 0.46, 1, 0.640, 0,1,0, 0.35400
0, 0.60, 9, 0.760, 1,0,0, 0.80800
1, 0.28, 0, 0.610, 0,0,1, 0.26800
0, 0.43, 5, 0.705, 0,1,0, 0.47600
1, 0.34, 4, 0.625, 1,0,0, 0.28300
0, 0.56, 6, 0.745, 0,0,1, 0.70400
1, 0.19, 2, 0.600, 0,1,0, 0.17400
0, 0.65, 8, 0.770, 1,0,0, 0.86800
1, 0.42, 3, 0.635, 0,0,1, 0.35900
0, 0.30, 7, 0.695, 0,1,0, 0.42200
1, 0.53, 1, 0.655, 1,0,0, 0.44800
0, 0.47, 9, 0.720, 0,0,1, 0.59400
1, 0.22, 0, 0.600, 0,1,0, 0.19800
0, 0.66, 5, 0.780, 1,0,0, 0.89200
1, 0.38, 4, 0.630, 0,0,1, 0.31800
0, 0.21, 6, 0.620, 0,1,0, 0.27100
1, 0.62, 2, 0.665, 1,0,0, 0.50300
0, 0.33, 8, 0.675, 0,0,1, 0.41200
1, 0.51, 3, 0.650, 0,1,0, 0.41500
0, 0.44, 7, 0.705, 1,0,0, 0.48800

Test data:

# people_test_40.txt
#
1, 0.27, 1, 0.620, 0,0,1, 0.28100
0, 0.58, 9, 0.750, 0,1,0, 0.75900
1, 0.32, 0, 0.625, 1,0,0, 0.25900
0, 0.24, 5, 0.640, 0,0,1, 0.32400
1, 0.60, 4, 0.660, 0,1,0, 0.52100
0, 0.39, 6, 0.685, 1,0,0, 0.42700
0, 0.26, 2, 0.690, 0,1,0, 0.40200
1, 0.49, 8, 0.645, 0,0,1, 0.38800
0, 0.53, 3, 0.730, 1,0,0, 0.61200
1, 0.30, 7, 0.615, 0,1,0, 0.23400
0, 0.61, 1, 0.760, 0,0,1, 0.78900
1, 0.22, 9, 0.610, 1,0,0, 0.21500
0, 0.45, 0, 0.700, 0,1,0, 0.49300
1, 0.57, 5, 0.670, 0,0,1, 0.51100
0, 0.28, 4, 0.695, 1,0,0, 0.40800
1, 0.39, 6, 0.630, 0,1,0, 0.31100
0, 0.51, 2, 0.720, 0,0,1, 0.57400
1, 0.33, 8, 0.625, 1,0,0, 0.26700
0, 0.64, 3, 0.770, 0,1,0, 0.83500
1, 0.18, 7, 0.600, 0,0,1, 0.16200
0, 0.42, 1, 0.695, 1,0,0, 0.46100
1, 0.59, 9, 0.680, 0,1,0, 0.56300
0, 0.36, 0, 0.685, 0,0,1, 0.43400
1, 0.23, 5, 0.610, 1,0,0, 0.22100
0, 0.67, 4, 0.780, 0,1,0, 0.92400
1, 0.45, 6, 0.645, 0,0,1, 0.37800
0, 0.31, 2, 0.700, 1,0,0, 0.43500
1, 0.52, 8, 0.655, 0,1,0, 0.43900
0, 0.20, 3, 0.620, 0,0,1, 0.26200
1, 0.48, 7, 0.650, 1,0,0, 0.37100
0, 0.55, 1, 0.740, 0,1,0, 0.63800
1, 0.29, 9, 0.620, 0,0,1, 0.27300
0, 0.63, 0, 0.765, 1,0,0, 0.80100
1, 0.41, 5, 0.640, 0,1,0, 0.34200
0, 0.22, 4, 0.630, 0,0,1, 0.28600
1, 0.64, 6, 0.675, 1,0,0, 0.61400
0, 0.37, 2, 0.690, 0,1,0, 0.44100
1, 0.50, 8, 0.650, 0,0,1, 0.42200
0, 0.27, 3, 0.695, 1,0,0, 0.39800
1, 0.44, 7, 0.640, 0,1,0, 0.33900

Data smoothing is a classical statistics technique that removes noise from a dataset so that overall patterns and trends are visible. There are many algorithms for smoothing. I put together a demo of data smoothing using the Nadaraya-Watson (NW) kernel smoothing technique, with the C# language.

The output of the demo program is:

Begin Nadaraya-Watson kernel smoothing demo

Loading simple XXy data
Done

First three data X:
  0.1000  0.1000
  0.2000  0.2000
  0.3000  0.3000

First three data y:
  0.0998
  0.3987
  0.0955

Creating NW kernel smoother with RBF sigma = 0.1600
Done

Computing smoothed y values
Done

   X             y        smooth y
-------------------------------------
  0.10  0.10    0.0998    0.2069
  0.20  0.20    0.3987    0.2303
  0.30  0.30    0.0955    0.2474
  0.50  0.50    0.5794    0.5599
  0.60  0.60    0.5646    0.6681
  0.70  0.70    0.9442    0.7357
  0.90  0.90    0.4833    0.5556
  1.00  1.00    0.5415    0.5255

End NW kernel smoothing demo

The source data is:

0.1000, 0.1000, 0.0998
0.2000, 0.2000, 0.3987
0.3000, 0.3000, 0.0955
0.5000, 0.5000, 0.5794
0.6000, 0.6000, 0.5646
0.7000, 0.7000, 0.9442
0.9000, 0.9000, 0.4833
1.0000, 1.0000, 0.5415

There are eight data items. The first two values on each line are the elements of a two-dimensional X vector, such as (0.10, 0.10). The third value is a dependent y value. The data is artificial but you can imagine the X values are inputs to some system and the y values are outputs.

Data smoothing is most often used with small datasets that have a single x value and a single y value. This allows the source data and the smoothed data to be visualized using a graph. However, I used data with two x values, such as (0.10, 0.10), where both values are artificially the same, to demonstrate that NW kernel smoothing can use multidimensional X data with several values.

Understanding the Radial Basis Function

In order to understand NW kernel smoothing, you must understand the radial basis function (RBF) function. The RBF function computes a measure of similarity between two vectors. There are several RBF versions. The version used by the demo is often called the sigma (bandwidth) version:

rbf(v1, v2) = exp(-1 * ||v1 - v2||^2 / (2 * sigma^2))

The v1 and v2 are vectors. The exp() function is Euler’s constant (approximately 2.71828) raised to a power. The ||v1 – v2||^2 is the squared Euclidean distance between vectors v1 and v2. The sigma, also called the bandwidth, is an arbitrary constant.

If the two vectors passed to RBF are the same, then rbf(v1, v2) = 1.0 (maximum similarity). As the two vectors get more and more different, the value of rbf(v1, v2) gets closer and closer to 0. RBF is symmetric so rbf(v1, v2) = rbf(v2, v1).

Suppose v1 = (2.50, -3.25, 1.00) and v2 = (2.0, -3.0, 1.0) and sigma = 1.5, then the squared distance term is:

||v1 - v2||^2 = (2.50 - 2)^2 + (-3.25 - (-3))^2 + (1.00 - 1)^2
              = (0.50)^2 + (-0.25)^2 + (0.00)^2
              = 0.2500 + 0.0625 + 0.0000
              = 0.3125

and the RBF value is:

rbf(v1, v2) = exp( -1 * ||v1 - v2||^2 / (2 * sigma^2) )
            = exp( -0.3125 / (2 * 1.5^2) )
            = exp( -0.3125 / 4.50 )
            = exp(-0.0694)

Computing a NW Kernel Smoothed y Value

Computing a NW smoothed data point is best explained using a concrete example. Suppose you want to compute the smoothed y’ value for x = (0.70, 0.70) in the demo data. The actual y value is 0.9442. The smoothed y’ value is 0.7357. It is computed as follows:

x = (0.70, 0.70) 

x0 = (0.10, 0.10)  w0 = rbf(x, x0) = 0.0000  y0 = 0.0998  w0 * y0 = 0.0000
x1 = (0.20, 0.20)  w1 = rbf(x, x1) = 0.0001  y1 = 0.3987  w1 * y1 = 0.0000
x2 = (0.30, 0.30)  w2 = rbf(x, x2) = 0.0019  y2 = 0.0955  w2 * y2 = 0.0001
x3 = (0.50, 0.50)  w3 = rbf(x, x3) = 0.2096  y3 = 0.5794  w3 * y3 = 0.1214
x4 = (0.60, 0.60)  w4 = rbf(x, x4) = 0.6766  y4 = 0.5646  w4 * y4 = 0.3820
x5 = (0.70, 0.70)  w5 = rbf(x, x5) = 1.0000  y5 = 0.9442  w5 * y5 = 0.9442
x6 = (0.90, 0.90)  w6 = rbf(x, x6) = 0.2096  y6 = 0.4833  w6 * y6 = 0.1013
x7 = (1.00, 1.00)  w7 = rbf(x, x7) = 0.0297  y7 = 0.5415  w7 * y7 = 0.0161
                                    -------                        -------
                       sum weights = 2.1276        sum weighted y = 1.5653

smooth y' = sum weighted y / sum weights
          = 1.5653 / 2.1276
          = 0.7357

There are eight data points. A weight, wi, is computed for each data point, where the weights are the RBF values between the x value and each data points. Each weight is multiplied by associated actual y value, and those products are summed to give the numerator term of the smoothed value. The denominator is the sum of the weight values.

I miss the days when I used to read a daily newspaper to get all my news and information. Here is some old health advice.

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

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

namespace KernelNWSmoothing
{
  internal class KernelNWSmoothingProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin Nadaraya-Watson" +
        " kernel smoothing demo ");

      Console.WriteLine("\nLoading simple XXy data ");
      string dataFile = "..\\..\\..\\Data\\simple_xxy_8.txt";
      double[][] X = MatLoad(dataFile,
        new int[] { 0, 1 }, ',', "#");
      double[] y = MatToVec(MatLoad(dataFile,
        new int[] { 2 }, ',', "#"));
      Console.WriteLine("Done ");

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

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

      double sigma = 0.16; // larger sigma  = smoother
      Console.WriteLine("\nCreating NW kernel smoother" +
        " with RBF sigma = " + sigma.ToString("F4"));
      KernelSmoother ks = new KernelSmoother("rbf", sigma);
      Console.WriteLine("Done ");

      Console.WriteLine("\nComputing smoothed y values ");
      ks.Fit(X, y); // copy dataX and dataY
      Console.WriteLine("Done ");

      Console.WriteLine("\n   X         y        smooth y");
      Console.WriteLine("--------------------------------");
      for (int i = 0; i "lt" X.Length; ++i)
      {
        double[] x = X[i];
        double actY = y[i];
        double smoothY = ks.SmoothedY(x);
        VecShow(x, 2, 6, false);
        Console.Write("    " + actY.ToString("F4"));
        Console.WriteLine("    " + smoothY.ToString("F4"));
      }

      Console.WriteLine("\nEnd NW kernel smoothing 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,
      bool newLine)
    {
      for (int i = 0; i "lt" vec.Length; ++i)
        Console.Write(vec[i].ToString("F" + dec).
          PadLeft(wid));
      if (newLine == true)
        Console.WriteLine("");
    }

  } // class Program

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

  public class KernelSmoother
  {
    public string kernel; // 
    public double sigma;  // RBF sigma, aka bandwidth
    public double[][] dataX;
    public double[] dataY;

    public KernelSmoother(string kernel, double sigma)
    {
      if (kernel != "rbf")
        Console.WriteLine("Only rbf kernel supported ");
      this.kernel = kernel; // dummy not used this version
      this.sigma = sigma;
    }

    public void Fit(double[][] dataX, double[] dataY)
    {
      this.dataX = dataX; // copy by reference
      this.dataY = dataY;
    }

    public double SmoothedY(double[] x)
    {
      int n = this.dataX.Length;
      double numer = 0.0; // sum of weighted y values
      double denom = 0.0; // sum of weights (RBF values)
      for (int i = 0; i "lt" n; ++i)
      {
        double currRbf = RbfSigma(x, this.dataX[i],
          this.sigma);
        numer += currRbf * this.dataY[i];
        denom += currRbf;
      }
      return numer / denom; // assume no div by 0
    }

    private static double RbfSigma(double[] x1, double[] x2,
      double sigma)
    {
      // rbf(x1,x2) = exp(-1 * ||x1-x2||^2 / (2 * sigma^2))
      int dim = x1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" dim; ++i)
        sum += (x1[i] - x2[i]) * (x1[i] - x2[i]);
      return Math.Exp(-1 * sum / (2 * sigma * sigma));
    }

    private static double RbfGamma(double[] x1, double[] x2,
      double gamma)
    {
      // rbf(x1,x2) = exp(-1 * gamma * ||x1-x2||^2
      int dim = x1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" dim; ++i)
        sum += (x1[i] - x2[i]) * (x1[i] - x2[i]);
      return Math.Exp(-1 * gamma * sum);
    }

  } // class KernelSmoother

} // ns

Demo data:

# simple_xxy_8.txt
#
0.1000, 0.1000, 0.0998
0.2000, 0.2000, 0.3987
0.3000, 0.3000, 0.0955
0.5000, 0.5000, 0.5794
0.6000, 0.6000, 0.5646
0.7000, 0.7000, 0.9442
0.9000, 0.9000, 0.4833
1.0000, 1.0000, 0.5415

Suppose you have a tiny dataset of nine items like this:

Sex Hand  Age  Height   City      Score (y)
-------------------------------------------
 M  right  24   tall    Anaheim    0.60		
 F  left   49   short   Anaheim    0.70
		
 M  right  33   medium  Boulder    0.20		
 M  ambi   38   short   Boulder    0.30		
 F  left   51   medium  Boulder    0.40
 M  right  28   tall    Boulder    0.50
		
 M  right  38   medium  Chicago    0.20		
 F  left   62   short   Chicago    0.30		
 M  ambi   27   tall    Chicago    0.40

Your goal is to predict Score from Sex, Hand, Age, Height, and City. And you want to use a tree-based regression system (bagging tree regression, random forest regression, AdaBoost.R2 regression, or gradient boosting regression).

Tree-Based Regression with Categorical Predictor Variables

The numeric Age variable can be left as-is, or it can be normalized (for example by dividing all age values by 100).

But the Sex, Hand, Height, and City values must be encoded into numeric values. (Note: Some library implementations of tree-based regression have built-in support for categorical predictors.)

The simplest possible approach is to always use ordinal encoding for any categorical predictor. Sometimes this works well, but sometimes a more sophisticated encoding approach works better.

Here is my standard set of five techniques to encode categorical variables for tree-based regression. Note that these are my personal heuristics, not some sort of generally accepted heuristics.

cardinality    encoding
================================================
natural order  plain ordinal
2 (binary)     0-1
3 to 4         one-hot
5 to 10        one-hot (or target-sorted ordinal)
10+            smoothed-target (aka impact)

The Sex variable is binary because it has exactly two possible values. Therefore it can be encoded as M = 0, F = 1 (or vice versa).

The Height variable has natural math order, with short is less than medium is less than tall. Therefore it can be plain ordinal encoded: short = 0, medium = 1, tall = 2.

The Hand variable has 3 possible values (where the values don’t have an inherent order) so it can be one-hot encoded: left = 1 0 0, right = 0 1 0, ambi (ambidextrous) = 0 0 1. (Where the order is arbitrary).

The City variable has just 3 possible values, so in a non-demo scenario it should be one-hot encoded. But suppose there were 5 through 10 possible cities. It could be encoded as target-sorted. The average of the target y values for each city is computed, the averages are sorted, and then the sort order is used as the encoding. For this data, the average y value for Anaheim is (0.60 + 0.70) / 2 = 0.65. For Boulder the average is 0.35. For Chicago the average is 0.30. If the averages are sorted from smallest to largest, the order is Chicago (0.30), then Boulder (0.35), then Anaheim (0.65). Therefore the target-sorted ordinal encoding is Chicago = 0, Boulder = 1, Anaheim = 2.

The target-sorted ordinal encoding is somewhat controversial. Many of my colleagues avoid it, and use one-hot encoding for up to 10 possible values.

Now suppose there were 10 or more cities. They would be smoothed-target encoded (also called impact encoded). The equation to encode a particular categorical value is:

x’ = [(n_cat * y_cat) + (m * y_global)] / (n_cat + m)

The n_cat is the number of times the value occurs in the dataset. The y_cat is the average of the y values for the value (just like target-sorting above). The y_global is the average of all the y values in the dataset. The m value is a hyperparameter. It is often set to 10, so I’ll use m = 10.

For this data, the overall average y value is 0.4000. The averages for Anaheim, Boulder, Chicago are 0.65, 0.35, 0.30 respectively. The counts of each values for Anaheim, Boulder, Chicago are 2, 4, 3 respectively.

Therefore, if m is set to 10, the encoding calculations are:

Anaheim: [(2 * 0.65) + (10 * 0.4000)] / (2 + 10) = 0.4417
Boulder: [(4 * 0.35) + (10 * 0.4000)] / (4 + 10) = 0.3857
Chicago: [(3 * 0.30) + (10 * 0.4000)] / (3 + 10) = 0.3769

Let me re-emphasize that I’m demonstrating smoothed-target (aka impact) encoding on just three values (rather than when there are 10 or more possible values) just to keep the example calculations simple.

In theory, the value of m should be determined using trial-and-error experimentation. But in practice, the value of 10 is often just used. Larger values of m produce stronger shrinkage toward the global mean, but the actual effect of using different values of m can vary depending on the exact nature of the dataset being used.

In practice, when using smoothed-target encoding, to prevent the model from memorizing rare values, you should consider grouping low-frequency categories (for example occurring less than x% of the time, where the x% depends on many factors) into an ‘Other’ bucket before encoding.

When using smoothed-target encoding, there is a danger of target leakage. Therefore, you should perform a train-test split first, use the training data to calculate the smoothed-target encoding, and then use those encoding values to transform the test data too. You should not encode the entire source dataset and then perform a train-test split.

In machine learning, data preparation is almost always by far the most time-consuming part of the process.

Most of my friends, and I, are fascinated by all kinds of strange things that don’t interest normal people. For reasons that are essentially unknown to me, I am intrigued by movies that feature a creature that makes noises that are translated to “chittering” in the closed caption text. The are an amazing number of chittering instances in science fiction movies, but you’ll only notice these chitterings if you’re actively looking for them.

Right: In “Project Hail Mary” (2026), the Sun is being eaten slowly but surely by some sort of space organisms. A crew of three travel to the distant star Tau Ceti to discover why that star is not being consumed. Only one member of the crew, Ryland, survives the trip. Upon arriving at Tau Ceti, Ryland discovers another spaceship from planet Erid with the same mission. Only one of the alien crew, Rocky, survived the trip. For the next 20 minutes or so, until Ryland and Rocky learn how to communicate, Rocky is in non-stop chittering mode. This movie got excellent reviews. I rarely agree with reviews but I really liked this movie when I watched it with my pal Ken. My personal grade = A.

Bottom line: Encoding categorical predictor variables for use in tree-based regression systems is a surprisingly complex topic. There are two main encoding schemes:

1. Always use simple ordinal encoding.
2. Use ordinal, or one-hot, or smoothed target encoding:
  a. binary predictor = use 0-1 encoding (ordinal encoding)
  b. categorical with natural ordering = use ordinal encoding
  c. categorical no ordering, 3-10 values = use one-hot encoding
  d. categorical no ordering, 10+ values = smooth target encoding

The first scheme is extremely simple. The second scheme is quite complicated. Sometimes the first, simple, encoding scheme works better. Sometimes the second, complex encoding scheme works better.

This blog post discusses the first, simple encoding scheme. The second, complex, encoding scheme is discussed in other blog posts. Search for “smoothed target”.

Suppose you want to predict a single numeric value (a regression problem) using a machine learning technique that is based on collections of decision trees — specifically, bagging regression, random forest regression, AdaBoost.R2 regression, or gradient boost regression. (The encoding techniques for math-based regression techniques such as kernel ridge regression and neural network regression is a completely different story).

And suppose you have a dataset of people that looks like:

Sex Age  State  Height  Income    Politics
-------------------------------------------------
M   41   GA     68.5   59100.00   moderate
F   23   AL     61.0   36800.00   conservative
F   18   CA     64.0   42600.00   liberal
M   56   CT     71.0   62400.00   moderate
. . .

Your goal is to predict Income from the other five variables. The numeric Age and Height predictors, and the target Income variable can be left as-is or normalized. But the categorical Sex, State, and Politics variables must be encoded. Encoding categorical variables for use with tree-based regression systems seems like it should be simple but in fact is an astonishingly complicated topic.

Note: There are some regression tree systems, such as CART trees, that have built-in support for categorical variables.

Example of Gradient Boost Regression on Ordinal Encoded Data

I will cut to the chase and say that, surprisingly, when all factors are taken into account, using simple ordinal encoding often works as well as the much more complicated ordinal, one-hot, smoothed target encoding. So for Sex, Male = 0, Female = 1. For State, Alabama = 0, Alaska = 1, . . Wyoming = 9. (Assuming only 9 States). For Politics, conservative = 0, moderate = 1, liberal = 2. Simple.

If numeric variables are divide-by-constant normalized, and the categorical variables are ordinal encoded, the data becomes:

Sex Age  State  Height  Income   Politics
-------------------------------------------------
0   0.41  9     0.685   0.59100   1
1   0.23  1     0.610   0.36800   0
1   0.18  4     0.640   0.42600   2
0   0.56  6     0.710   0.62400   1
. . .

Ordinal encoding imposes an arbitrary math ordering such as, “Wyoming is greater than California”. But a tree splits data using simple threshold checks (e.g., State less-than-or-equal 14). It can chain these splits down a branch to isolate specific numbers or groups of numbers, effectively carving out individual states regardless of their arbitrary integer values.

A simple, but less obvious alternative is count-encoding. You encode according to the counts of each value. For example, if the dataset has 324 items/people from California (because it’s a populous state) but just 9 people from Wyoming, then California = 324 and Wyoming = 9.

Count encoding can be better than ordinal encoding when the frequency of occurrence contains meaningful real-world signal that directly correlates with your target variable, in this case State population size, which might be a predictor of Income.

A huge advantage of ordinal-encoding and count-encoding is that you can encode categorical variables in the source dataset and then split into training and test sets, as opposed to many encoding techniques such as “target-mean with OOF and smoothing”, where you must split first, then fold, then encode — a real nightmare in practice. Put another way, there is no data leakage with ordinal encoding and count-encoding.

The target-mean encoding with OOF and smoothing technique often appeals to beginners because it looks very impressive with all kinds of equations and related math concepts. But in a production environment, sometimes simplicity trumps fancy math.

Machine learning encoding transforms categorical data into numeric data.

Right: Stone Boy could . . . well, turn himself into stone. But then he couldn’t move or speak. He applied to the Legion in Adventure Comics #306 (March 1963) but was rejected. Other hilarious rejectees included Antennae Boy (radio antennae on his head and radio in his brain), Chlorophyll Kid (make plants grow fast), Camera Kid (a video camera in his forehead), Flat Boy (make himself flat like a rug), Leadfeather Lad (make bird feathers very heavy), and Spider Girl (expanding hair). In a later issue, Adventure Comics #315 (December 1963), Stone Boy showed great courage and was offered membership in the Legion, but he declined so he could remain a member of the Legion of Substitute Heroes.

As usual with complex problems, explaining what the problem is, is more difficult than explaining the solution. So bear with me. But the bottom line is that I tried an experiment that worked, but gave only a small 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 Gradient Boosting regression technique uses a collection of simple decision trees — they’re called the learners or the estimators. Each tree is constructed sequentially, based on error information of the previous tree. Each tree gets slightly better. The final prediction is based on the sum of indirect predictions of the trees (it’s a bit complicated to explain).

Although Gradient Boosting 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 another long and complicated story.

I decided to implement Gradient Boost regression using Blind Trees learners. Bottom line: For my demo dataset, the technique worked surprisingly well.

Blind Trees are even weaker than standard decision trees. When constructing a standard decision tree, for each tree node, an optimal split column and split value are found. When constructing a Blind Trees tree, for each node, a split value is selected blindly, without looking at the target y values, and without computing variance reduction. This makes blind Trees very fast to construct.

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 C# Gradient Boosting with Blind Trees demo

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

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

First three train y:
  0.4840
  0.1568
  0.8054

Setting numTrees = 600
Setting maxDepth = 6
Setting minSamples = 2
Setting minLeaf = 1
Setting lrnRate = 0.0500

Creating and training GradientBoostRegression model
Done

Evaluating model

Accuracy train (within 0.10) = 1.0000
Accuracy test (within 0.10) = 0.8250

MSE train = 0.0000
MSE test = 0.0013

R2 train = 1.0000
R2 test = 0.9533

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

Initial prediction: 0.3493
t =   0  pred_res =  0.1700  delta =  0.0085  pred =  0.3578
t =   1  pred_res =  0.1785  delta =  0.0089  pred =  0.3668
t =   2  pred_res =  0.1230  delta =  0.0061  pred =  0.3729
t =   3  pred_res =  0.2392  delta =  0.0120  pred =  0.3849
t =   4  pred_res =  0.2090  delta =  0.0105  pred =  0.3953
. . .
t = 598  pred_res = -0.0001  delta = -0.0000  pred =  0.4844
t = 599  pred_res = -0.0001  delta = -0.0000  pred =  0.4844
Predicted y = 0.4844

End demo

A very interesting investigation.

Gradient Boost regression code is a wrapper around a collection of weak learners. I’m a big fan of science fiction movies of the 1950s. Several of these movies had aliens that were obviously just actors wrapped in costumes of some kind, but this didn’t reduce the appeal for me. Here are two of my favorites.

Right: In “This Island Earth” (1955), human-like aliens from Metaluna come to Earth to seek scientific help in their war against the planet of Zagon. The Metaluna aliens use 7-foot tall mutants as general purpose helpers and for security. My grade = A-.

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 GradientBoostWithBlindTrees
{
  internal class GradientBoostWithBlindTreesProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# Gradient Boosting" +
        " with Blind Trees demo ");

      // 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 };
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        [5], ',', "#"));

      string testFile =
        "..\\..\\..\\Data\\synthetic_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        [5], ',', "#"));
      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));

      int numTrees = 600;
      int maxDepth = 6;
      int minSamples = 2;
      int minLeaf = 1;
      int numSplitCols = -1; // use all columns
      double lrnRate = 0.05;  // 0.9900 .8250  0  .0013 

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

      Console.WriteLine("\nCreating and training" +
        " GradientBoostRegression model ");
      GradientBoostRegressor gbr =
        new GradientBoostRegressor(numTrees, maxDepth,
        minSamples, minLeaf, numSplitCols, lrnRate, 0);
      gbr.Train(trainX, trainY);
      Console.WriteLine("Done ");

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

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

      double r2Train = gbr.R2(trainX, trainY);
      Console.WriteLine("\nR2 train = " +
        r2Train.ToString("F4"));
      double r2Test = gbr.R2(testX, testY);
      Console.WriteLine("R2 test = " +
        r2Test.ToString("F4"));

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

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

    // ------------------------------------------------------
    // helpers for Main()
    // ------------------------------------------------------

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

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

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

  } // class Program

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

  public class GradientBoostRegressor
  {
    public double lrnRate;
    public int nTrees;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    public int numSplitCols;
    public List"lt"BlindTreeRegressor"gt" trees;
    public double pred0;  // initial prediction
    private Random rnd;

    public GradientBoostRegressor(int nTrees, int maxDepth,
      int minSamples, int minLeaf, int numSplitCols,
      double lrnRate, int seed=0)
    {
      this.nTrees = nTrees;
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;
      this.trees = new List"lt"BlindTreeRegressor"gt"();
      this.lrnRate = lrnRate;
      this.rnd = new Random(seed);
    }

    public void Train(double[][] trainX, double[] trainY)
    {
      int n = trainX.Length;
      this.pred0 = Mean(trainY);

      double[] preds = new double[n]; //each data item
      for (int i = 0; i "lt" n; ++i)
        preds[i] = this.pred0;

      for (int t = 0; t "lt" this.nTrees; ++t) // each tree
      {
        double[] residuals = new double[n]; // for curr tree
        for (int i = 0; i "lt" n; ++i)
          residuals[i] = trainY[i] - preds[i];

        BlindTreeRegressor bt =
          new BlindTreeRegressor(this.maxDepth,
          this.minSamples, this.minLeaf, this.numSplitCols,
          false, this.rnd.Next(0,1_000_000));
        bt.Train(trainX, residuals); // predict residuals

        for (int i = 0; i "lt" n; ++i)
        {
          double predResidual = bt.Predict(trainX[i]);
          preds[i] += this.lrnRate * predResidual;
        }
        this.trees.Add(bt);
      }
    } // Train

    public double Predict(double[] x, bool verbose = false)
    {
      double result = this.pred0;
      if (verbose == true)
      {
        Console.WriteLine("\nInitial prediction: " +
        result.ToString("F4"));
      }
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;

        if (verbose == true)
        {
          if (t "gte" 0 && t "lte" 4 ||
            t "gte" this.nTrees - 2 &&
            t "lte" this.nTrees - 1)
          {
            Console.Write("t = " + t.ToString().PadLeft(3) +
              "  pred_res = " + predResidual.ToString("F4").
              PadLeft(7));
            Console.Write("  delta = " + delta.ToString("F4").
              PadLeft(7));
            Console.WriteLine("  pred = " +
              result.ToString("F4").PadLeft(7));
          }
          if (t == 5)
            Console.WriteLine(". . . ");
        }
      }
      return result;
    }

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

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

    public 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)
        sum += dataY[i];
      double meanActual = sum / n;

      double ssRes = 0.0;
      double ssTot = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes += (dataY[i] - predY) * (dataY[i] - predY);
        ssTot += (dataY[i] - meanActual) *
          (dataY[i] - meanActual);
      }
      if (Math.Abs(ssTot) "lt" 1.0e-12)
        return 0.0;
      else
        return 1.0 - (ssRes / ssTot);
    }

    private static double Mean(double[] data)
    {
      int n = data.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += data[i];
      return sum / n;
    }

  } // class GradientBoostRegression

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

  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;   // Explicitly isolate node boundaries
          currNode.right = -1;
          currNode.colIdx = -1;
          continue;
        }

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

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

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

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

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

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

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

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

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

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

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

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

        if (currNode == null) break;

        lastValidValue = currNode.value;

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

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

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

    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 “Neural Network Bistratum Regression Using C#” in the October 2026 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2026/10/01/neural-network-bistratum-regression-using-csharp.aspx.

Neural network bistratum regression is a machine learning technique that uses a neural network with exactly two hidden layers to predict a single numeric value.

A simple neural network has one hidden layer of processing nodes. The term deep neural network means a neural network that has two or more hidden layers. In practice, the term deep neural network is usually applied to a network used for image processing or natural language processing, and which has hundreds of layers. For convenience, I use the non-standard term bistratum (“two layers”, borrowed from biology) to indicate a neural network that has 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 five 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.

Neural network weights and biases are computed by looking at training data with known input values and known correct output values. The idea is to find values of the weights and biases so that computed output values are as close as possible to the target correct values. Put another way, the values of the weights and biases are those that minimize the mean squared error (or some other loss metric) between computed and target values.

There are several algorithms to find the values of neural network weights and biases. For relatively simple neural networks, the most common technique is called stochastic gradient descent (SGD). There are many advanced techniques that are related to SGD including Adam (adaptive moment estimation), RMSProp (root mean squared propagation), and others.

The demo program can be used as a template for most regression problems. The architecture parameter to explore is the number of nodes in the two hidden layers. The training hyperparameters to explore are the max epochs, learning rate, and decay.

I have written software for over 50 years. My experience gives me a pretty good subjective notion of what makes a beautiful machine learning model, and what makes an ugly model.

I better stick to evaluating software, not fashion.

The form of a linear regression prediction equation is y’ = (w0 * x0) + (w1 * x1) + . . + (wn * xn) + b where y’ is the predicted value, the xi are predictor values, the wi are constants called model weights, and b is a constant called the bias. For example, y’ = predicted bank account balance = (-0.54 * age) + (0.38 * height) + (0.11 * experience) + 0.72.

Training the model is the process of finding the values of the weights and bias so that predicted y values are close to known correct target y values in a set of training data.

There are several ways to train a linear regression model. Four of the most common techniques are iterative stochastic gradient descent (SGD) training, closed-form training using explicit matrix pseudo-inverse, closed-form training using implicit matrix pseudo-inverse via a back-substitution solve function, and closed-form training using explicit matrix left pseudo-inverse.

One day, I decided to implement the left pseudo-inverse training technique. Briefly, the math equation for explicit pseudo-inverse training is:

w = pinv(X) * y

Here w is the vector of the weights and the bias you are trying to find. The pinv() function is the Moore-Penrose pseudo-inverse of a matrix. X is a design matrix of the training X predictors of the training data. The * is matrix-to-vector multiplication, and y is a vector of target y values of the training data. The pinv() function can be implemented using singular value decomposition (SVD) or QR decomposition.

A design matrix adds a leading column of 1s, which takes the bias into account. For example, if a set of training predictors stored in a matrix X is:

0.50  0.77  0.32
0.92  0.41  0.84
0.64  0.73  0.43

The associated design matrix is:

1.0  0.50  0.77  0.32
1.0  0.92  0.41  0.84
1.0  0.64  0.73  0.43

By explicit pseudo-inverse training, I mean you use an actual pinv() Moore-Penrose pseudo-inverse function. But instead of using an explicit pseudo-inverse function, it’s possible to simulate pseudo-inverse using a regular matrix inverse function. The math equation of the left pseudo-inverse form of training is:

w = inv(Xt * X) * Xt * y

Here inv() is a any matrix inverse function. Usually Cholesky inverse is used, but any general purpose algorithm such as LUP inverse, or QR inverse, or SVD inverse, or Newton inverse can be used. X is the design matrix, Xt is the transpose of X (rows and columns exchanged). This left pseudo-inverse approach is simpler than using explicit pseudo-inverse, but X * Xt must be invertible, which you won’t know in advance. One way of dealing with this issue is to “condition” the Xt * X matrix by adding a small constant (about 1.0e-6) to the diagonal elements. The Xt * X matrix will be positive definite, so instead of using a general purpose matrix inverse technique (of which there are about a dozen techniques), you can use a simpler matrix inverse based on Cholesky decomposition.

Before I go any further, let me point out that the left pseudo-inverse technique only works for small X matrices. The problem is not the Cholesky inverse, the problem is computing Xt * X. If X has n rows, there will be roughly n^2 multiplications and additions. This could be millions of operations and arithmetic overflow or underflow could easily happen.

I put together a demo. For the demo, I used one of my standard sets of synthetic data that 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 five values on each line are the predictors (sometimes called features). The last value on each line is the target y value to predict. There are 200 training items and 40 test items.

The output of my demo is:

Begin C# linear regression using left
  pseudo-inverse training

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

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

First three train y:
  0.4840
  0.1568
  0.8054

Creating and training Linear Regression model
  using left p-inverse
Done

Coefficients/weights:
-0.2656  0.0333  -0.0454  0.0358  -0.1146
Bias/constant: 0.3619

Evaluating model

Accuracy train (within 0.10) = 0.4600
Accuracy test (within 0.10) = 0.6500

MSE train = 0.0026
MSE test = 0.0020

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

Predicted y = 0.5329

End demo

One of the minor design problems when using left pseudo-inverse training is that you have to decide which one of the dozens of possible matrix inverse algorithms to use. For my demo, instead of using Cholesky inverse, just for hoots I arbitrarily decided to use QR decomposition, Householder algorithm. The key statements in the Train() method are:

double[][] X = MatToDesign(trainX);
double[][] Xt = MatTranspose(X);
double[][] XtX = MatProduct(Xt,X);  // this could fail
// optionally condition Xt*X here
double[][] invXtX = MatInverseQR(XtX);  // this could fail
double[][] invXtXXt = MatProduct(invXtX, Xt);
double[] biasAndWts = MatVecProd(invXtXXt, trainY);

The corresponding statements in an explicit pseudo-inverse Train() method would be

double[][] X = MatToDesign(trainX);  // design X
double[][] Xpinv = MatPinv(X); // explicit pseudo-inverse
double[] biasAndWts = MatVecProd(Xpinv, trainY);

An interesting exploration. Again: the code in this blog post is not for real use in a production environment.

Using implicit pseudo-inverse training or explicit pseudo-inverse training is to some extent a matter of style.

I know very little about art, and I know even less about Japanese art. I went to Tokyo a few month ago. Before going, I read a bit about things Japanese. Japanese bijinga is an art genre that translates to “pictures of beautiful women”. Within that genre, different artists have very different styles. Here are three examples of modern bijinga art. Left: By artist Ayana Otake. Center: By artist Yasunari Ikenaga. Right: By artist Ichiro Tsuruta.

Demo code. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols.

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

namespace LinearRegressionPseudoInverseImplicitQR
{
  internal class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression" +
        " using implicit pseudo-inverse (QR) training ");

      // 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 };
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        new int[] { 5 }, ',', "#"));

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

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

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

      // 2. create and train implicit pseudo-inverse
      Console.WriteLine("\nCreating and training" +
        " Linear Regression model using implied p-inverse ");
      LinearRegressor model = new LinearRegressor();
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");

      // 2b. show model parameters
      Console.WriteLine("\nCoefficients/weights: ");
      for (int i = 0; i "lt" model.weights.Length; ++i)
        Console.Write(model.weights[i].ToString("F4") + "  ");
      Console.WriteLine("\nBias/constant: " +
        model.bias.ToString("F4"));

      // 3. evaluate model
      Console.WriteLine("\nEvaluating 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 predict first training item
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 9);
      double predY = model.Predict(x);
      Console.WriteLine("\nPredicted 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[][] mat)
    {
      int nRows = mat.Length;
      int nCols = mat[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = mat[i][j];
      return result;
    }

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

  } // class Program

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

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

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

    public LinearRegressor(int seed = 0)  // ctor
    {
      this.weights = new double[0];
      this.bias = 0;
      this.rnd = new Random(seed); // no need this version
    }

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

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

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

    public void Train(double[][] trainX, double[] trainY)
    {
      // implicit pseudo-inverse training
      // w = inv(Xt * X) * Xt * y
      int dim = trainX[0].Length;
      this.weights = new double[dim];

      double[][] X = MatToDesign(trainX);
      double[][] Xt = MatTranspose(X);
      double[][] XtX = MatProduct(Xt, X);
      for (int i = 0; i "lt" XtX.Length; ++i)
        XtX[i][i] += 1.0e-6;  // condition Xt*X
      double[][] invXtX = MatInverseQR(XtX);
      double[][] invXtXXt = MatProduct(invXtX, Xt);
      double[] biasAndWts = MatVecProd(invXtXXt, trainY);
      this.bias = biasAndWts[0];
      for (int i = 1; i "lt" biasAndWts.Length; ++i)
        this.weights[i - 1] = biasAndWts[i];
      return;
    } // Train()

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

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

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

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

    // ------------------------------------------------------
    // helpers
    // ------------------------------------------------------

    private static double[][] MatInverseQR(double[][] M)
    {
      // if A = Q * R, Ainv = Rinv * Qinv
      double[][] Q; double[][] R;
      MatDecomposeQR(M, out Q, out R, reduced: true);

      double[][] Rinv = MatInvUpperTri(R);
      double[][] Qinv = MatTranspose(Q);
      double[][] result = MatProduct(Rinv, Qinv);
      return result;
    }

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

    private static void MatDecomposeQR(double[][] M,
      out double[][] Q, out double[][] R, bool reduced)
    {
      // QR decomposition, Householder algorithm.
      // see rosettacode.org/wiki/QR_decomposition
      int m = M.Length;
      int n = M[0].Length;

      if (m "lt" n)
        throw new Exception("No rows less than cols");

      double[][] QQ = MatIdentity(m); // working Q
      double[][] RR = MatCopy(M); // working R

      int end;
      if (m == n) end = n - 1;
      else end = n;

      for (int i = 0; i "lt" end; ++i)
      {
        double[][] H = MatIdentity(m);
        double[] a = new double[m - i]; // corr
        int k = 0;
        for (int ii = i; ii "lt" m; ++ii) // corr
          a[k++] = RR[ii][i];

        double normA = VecNorm(a);
        if (a[0] "lt" 0.0 && normA "gt" 0.0) // corr
          normA = -normA;
        else if (a[0] "gt" 0.0 && normA "lt" 0.0)
          normA = -normA;

        double[] v = new double[a.Length];
        for (int j = 0; j "lt" v.Length; ++j)
          v[j] = a[j] / (a[0] + normA);
        v[0] = 1.0;

        // Householder algorithm
        double[][] h = MatIdentity(a.Length);
        double vvDot = VecDot(v, v);
        double[][] A = VecToMat(v, v.Length, 1);
        double[][] B = VecToMat(v, 1, v.Length);
        double[][] AB = MatProduct(A, B);

        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            h[ii][jj] -= (2.0 / vvDot) * AB[ii][jj];

        // copy h[][] into lower right corner of H[][]
        int d = m - h.Length; // corr
        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            H[ii + d][jj + d] = h[ii][jj];

        QQ = MatProduct(QQ, H);
        RR = MatProduct(H, RR);
      } // i

      if (reduced == false)
      {
        Q = QQ; // working results into the out params
        R = RR;
        return;
      }
      //else if (reduced == true)
      {
        int qRows = QQ.Length; int qCols = QQ[0].Length;
        int rRows = RR.Length; int rCols = RR[0].Length;
        // assumes m "gte" n !!

        // square-up R
        int dim = Math.Min(rRows, rCols);
        double[][] Rsquared = MatMake(dim, dim);
        for (int i = 0; i "lt" dim; ++i)
          for (int j = 0; j "lt" dim; ++j)
            Rsquared[i][j] = RR[i][j];

        // Q needs same number columns as R
        // so that inv(R) * trans(Q) works
        double[][] Qtrimmed = MatMake(qRows, dim);
        for (int i = 0; i "lt" qRows; ++i)
          for (int j = 0; j "lt" dim; ++j)
            Qtrimmed[i][j] = QQ[i][j];

        Q = Qtrimmed;
        R = Rsquared;
        return;
      }
    } // MatDecomposeQR()

    // ------------------------------------------------------
    
    private static double[] MatVecProd(double[][] M,
      double[] v)
    {
      // return a regular vector
      int nRows = M.Length;
      int nCols = M[0].Length;
      int n = v.Length;
      if (nCols != n)
        throw new Exception("non-comform in MatVecProd");

      double[] result = new double[nRows];
      for (int i = 0; i "lt" nRows; ++i)
        for (int k = 0; k "lt" nCols; ++k)
          result[i] += M[i][k] * v[k];

      return result;
    }

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

    private static double[][] MatToDesign(double[][] M)
    {
      // add a column of 1s
      int nRows = M.Length;
      int nCols = M[0].Length;
      double[][] result = new double[M.Length][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols + 1];

      for (int i = 0; i "lt" nRows; ++i)
      {
        result[i][0] = 1.0;
        for (int j = 1; j "lt" nCols + 1; ++j)
          result[i][j] = M[i][j - 1];
      }

      return result;
    }

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

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

    private static double[][] MatInvUpperTri(double[][] U)
    {
      int n = U.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[j][k] -= result[j][i] * U[i][k];
          }
          result[j][k] /= (U[k][k] + 1.0e-8); // avoid 0
        }
      }
      return result;
    }

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

    private static double[][] MatIdentity(int n)
    {
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      return result;
    }

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

    private static double[][] MatProduct(double[][] A,
      double[][] B)
    {
      int aRows = A.Length;
      int aCols = A[0].Length;
      int bRows = B.Length;
      int bCols = B[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

      double[][] result = MatMake(aRows, bCols);

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

    private static double[][] MatTranspose(double[][] M)
    {
      int nr = M.Length;
      int nc = M[0].Length;
      double[][] result = MatMake(nc, nr);  // note
      for (int i = 0; i "lt" nr; ++i)
        for (int j = 0; j "lt" nc; ++j)
          result[j][i] = M[i][j];
      return result;
    }

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

    private static double[][] MatCopy(double[][] M)
    {
      int nr = M.Length; int nc = M[0].Length;
      double[][] result = MatMake(nr, nc);
      for (int i = 0; i "lt" nr; ++i)
        for (int j = 0; j "lt" nc; ++j)
          result[i][j] = M[i][j];
      return result;
    }

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

    private static double VecNorm(double[] vec)
    {
      int n = vec.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += vec[i] * vec[i];
      return Math.Sqrt(sum);
    }

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

    private static double[][] VecToMat(double[] vec,
        int nRows, int nCols)
    {
      double[][] result = MatMake(nRows, nCols);
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = vec[k++];
      return result;
    }

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

    private static double VecDot(double[] v1, double[] v2)
    {
      double result = 0.0;
      int n = v1.Length;
      for (int i = 0; i "lt" n; ++i)
        result += v1[i] * v2[i];
      return result;
    }

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

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