James D. McCaffrey

Logistic regression is a technique for binary classification. For example, you might want to predict the sex of a person (0 = male, 1 = female) based on age, State, income, and political leaning. I enjoy writing code, so one rainy evening (I live in the Pacific Northwest so this wasn’t unusual), while I was waiting for the rain to stop so I could walk my dogs, I decided to implement logistic regression, from scratch, using C#.

In logistic regression, the prediction output is a pseudo-probability value between 0.0 and 1.0. If the p-value is less than 0.5, the prediction is class 0. If the p-value is greater than or equal to 0.5, the prediction is class 1. For any problem, the choice of which label is 0 and which is 1, is arbitrary.

The prediction equation for logistic regression is:

z = (wo * x0) + (w1 * x1) + . . . + b
p = 1.0 / (1.0 + exp(-z))

The xi are predictor values. The wi are model weights. The b is a bias constant. Training the model is the process of finding values for the weights and the bias so that computed output probabilities closely match the target y (0 / 1) values.

For my demo, I used one of my standard artificial datasets. The raw data looks like:

F, 24, michigan, 29500.00, liberal
M, 39, oklahoma, 51200.00, moderate
F, 63, nebraska, 75800.00, conservative
. . .

The normalized and encoded data looks like:

1, 0.24, 1,0,0, 0.2950, 0,0,1
0, 0.39, 0,0,1, 0.5120, 0,1,0
1, 0.63, 0,1,0, 0.7580, 1,0,0
. . .

There are only three State values. The output of my demo is:

Begin C# linear regression demo with Newton-Raphson training

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

First three train X:
  0.2400  1.00  0.00  0.00  0.2950  0.00  0.00  1.00
  0.3900  0.00  0.00  1.00  0.5120  0.00  1.00  0.00
  0.6300  0.00  1.00  0.00  0.7580  1.00  0.00  0.00

First three train y:
  1.0000
  0.0000
  1.0000

Creating Logistic Regression model
Done

Start Newton-Raphson train
epoch =     0  MSE =   0.1274
epoch =    40  MSE =   0.1145
epoch =    80  MSE =   0.1145
epoch =   120  MSE =   0.1145
epoch =   160  MSE =   0.1145
Done

Weights/coefficients:
81.4712 1.8590 1.8903 0.2694 -80.7779 3.9275 1.1614 -1.0702
Bias/constant: 4.0187

Evaluating model

Accuracy train = 0.8450
Accuracy test = 0.7250

MSE train = 0.1145
MSE test = 0.1305

Confusion matrix:
 actual   0 |     17     9  | 0.6538
 actual   1 |      2    12  | 0.8571
----------
predicted:         0     1  | 0.7250

End demo

There are dozens of design alternatives for logistic regression. I used Newton-Raphson training with binary cross entropy (BCE) loss gradient updates, combined with mean squared error loss progress messages.

Newton-Raphson training is appealing because it does not require any training parameters such as a learning rate and maximum epochs, that must be determined by trial and error. The main alternative, stochastic gradient descent (SGD), does require a learning rate and max epochs. Newton-Raphson computes an nxn Hessian matrix, which isn’t feasible for extrememly large datasets, but Newton-Raphson is very fast compared to SGD training. L-BFGS training is possible too, but it isn’t well-suited for scratch implementation.

I implemented code to compute and display a confusion matrix. This is more or less standard procedure for binary classification problems.

I validated my logistic regression system by running the demo data through the scikit-learn library LogisticRegression module. I got essentially the same results. See https://jamesmccaffreyblog.com/2023/02/10/logistic-regression-using-the-scikit-library-in-visual-studio-magazine/.

My implementation can be improved. Currently, the training loop iterates a fixed number of iterations (maxEpochs). It’s possible to monitor the change in model weights and bias, and early-exit training when the largest change is less than some small value, like 1.0e-6 or so. But that’s a task for another day.

Good fun on a rainy Pacific Northwest night.

I am a big fan of science fiction films of the 1950s. Some of these movies have really bad special effects, but pretty good acting and story, make them entertaining (to me anyway).

Right: In “The Giant Claw” (1957), a giant bird creature somehow manages to get to Earth from an antimatter galaxy, which makes it invulnerable to all missiles and weapons. Scientists devise a way to remove the creature’s antimatter protection, and kill it with rockets. This is a pretty good movie but the production ran out of money before the special effects could be added, and resorted to low-budget Mexican special effects company. My grade for the movie = B.

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 LogisticRegressionPeopleNR
{
  internal class LogisticRegressionNRProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression" +
        " demo with Newton-Raphson training ");
      // sex age State income politics
      // 1, 0.24, 1,0,0, 0.2950, 0,0,1
      // 0, 0.39, 0,0,1, 0.5120, 0,1,0
      // 0    1   2 3 4     5    6 7 8
      //
      // predict sex from age State income politics

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

      string testFile =
        "..\\..\\..\\Data\\people_test_40.txt";
      double[][] testX =
        MatLoad(testFile, colsX, ',', "#");
      double[] testY =
        MatToVec(MatLoad(testFile,
        new int[] { 0 }, ',', "#"));
      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 model
      Console.WriteLine("\nCreating Logistic " +
        "Regression model ");
      LogisticRegression model =
        new LogisticRegression();
      Console.WriteLine("Done ");

      Console.WriteLine("\nStart Newton-Raphson train ");
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");

      // 2b.show model parameters
      Console.WriteLine("\nWeights/coefficients: ");
      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);
      Console.WriteLine("\nAccuracy train = " +
        accTrain.ToString("F4"));
      double accTest = model.Accuracy(testX, testY);
      Console.WriteLine("Accuracy test = " +
        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"));

      Console.WriteLine("\nConfusion matrix: ");
      int[][] cm = model.ConfusionMatrix(testX, testY);
      model.ShowConfusion(cm);

      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[][] M)
    {
      int nRows = M.Length;
      int nCols = M[0].Length;
      double[] result = new double[nRows * nCols];
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[k++] = M[i][j];
      return result;
    }

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

  } // class Program

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

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

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

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

    public int Predict(double[] x)
    {
      double p = ComputeP(x);
      if (p "lt" 0.5)
        return 0;
      else
        return 1;
    }

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

    public double ComputeP(double[] x)
    {
      double z = 0.0;
      for (int i = 0; i "lt" x.Length; ++i)
        z += this.weights[i] * x[i];
      z += this.bias;
      double p = this.Sigmoid(z);
      return p;
    }

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

    public double Sigmoid(double z)
    {
      if (z "lt" -6.0)
        return 0.0;
      else if (z "gt" 6.0)
        return 1.0;
      else
        return 1.0 / (1.0 + Math.Exp(-z));
    }

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

    public void Train(double[][] trainX, double[] trainY)
    {
      // Newton-Raphson technique
      // initialize weights and bias
      int dim = trainX[0].Length;
      int n = trainX.Length;
      int m = dim + 1;  // include bias as final parameter

      this.weights = new double[dim];
      this.bias = 0.0;
      int maxEpochs = trainX.Length; // heuristic

      // Newton-Raphson iterations
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        // gradient and Hessian
        double[] grad = new double[m]; // features + bias
        double[][] H = new double[m][]; // H is mxm
        for (int i = 0; i "lt" m; ++i)
          H[i] = new double[m];

        // accumulate gradient and Hessian
        for (int i = 0; i "lt" n; ++i) // each train item
        {
          double[] x = trainX[i];
          double y = trainY[i];

          double p = this.ComputeP(x);
          double error = p - y; // order is important
          double v = p * (1.0 - p);

          // include bias as a final feature
          double[] xBias = new double[m];
          for (int j = 0; j "lt" dim; ++j)
            xBias[j] = x[j];
          xBias[dim] = 1.0;

          // gradient
          for (int j = 0; j "lt" m; ++j)
            grad[j] += error * xBias[j];

          // Hessian
          for (int j = 0; j "lt" m; ++j)
            for (int k = 0; k "lt" m; ++k)
              H[j][k] += v * xBias[j] * xBias[k];
        }

        // average gradient and Hessian
        for (int j = 0; j "lt" m; ++j)
        {
          grad[j] /= n;
          for (int k = 0; k "lt" m; ++k)
            H[j][k] /= n;
        }

        // small diagonal for numerical stability
        for (int j = 0; j "lt" m; ++j)
          H[j][j] += 1.0e-8;

        // Cholesky decomposition: H = L * L^T
        double[][] L = new double[m][];
        for (int i = 0; i "lt" m; ++i)
          L[i] = new double[m];

        for (int i = 0; i "lt" m; ++i)
        {
          for (int j = 0; j "lte" i; ++j)
          {
            double sum = H[i][j];

            for (int k = 0; k "lt" j; ++k)
              sum -= L[i][k] * L[j][k];

            if (i == j)
              L[i][j] = Math.Sqrt(sum);
            else
              L[i][j] = sum / L[j][j]; // not 0
          }
        } // each item

        // solve L * z = -gradient
        // (forward substitution)
        double[] z = new double[m];

        for (int i = 0; i "lt" m; ++i)
        {
          double sum = -grad[i];
          for (int j = 0; j "lt" i; ++j)
            sum -= L[i][j] * z[j];
          z[i] = sum / L[i][i];
        }

        // Solve L^T * delta = z
        // (back substitution)
        double[] delta = new double[m];

        for (int i = m - 1; i "gte" 0; --i)
        {
          double sum = z[i];
          for (int j = i + 1; j "lt" m; ++j)
            sum -= L[j][i] * delta[j];
          delta[i] = sum / L[i][i];
        }

        // update weights and bias
        for (int j = 0; j "lt" dim; ++j)
          this.weights[j] += delta[j];
        this.bias += delta[dim];

        // show progress
        if (epoch % Math.Max(1, maxEpochs / 5) == 0)
        {
          double mse = this.MSE(trainX, trainY);
          Console.WriteLine("epoch = " +
            epoch.ToString().PadLeft(5) +
            "  MSE = " + mse.ToString("F4").PadLeft(8));
        }
      } // each epoch
    } // Train()

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

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

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

    public double Accuracy(double[][] dataX, double[] dataY)
    {
      int n = dataX.Length;
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" n; ++i)
      {
        double[] x = dataX[i];
        double y = (int)(dataY[i]);  // 0 or 1
        int predY = this.Predict(x);  // 0 or 1

        if (predY == y)
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numCorrect + numWrong);
    }

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

    public int[][] ConfusionMatrix(double[][] dataX,
      double[] dataY)
    {
      // assumes binary classification
      int[][] result = new int[2][];  // 2x2
      result[0] = new int[2];
      result[1] = new int[2];
      int n = dataX.Length;
      for (int i = 0; i "lt" n; ++i) // each data item
      {
        int predY = this.Predict(dataX[i]); // 0 or 1
        int actualY = (int)(dataY[i]);  // 0 or 1

        ++result[actualY][predY];

        //if (actualY == 0 && predY == 0)
        //  ++result[0][0];
        //else if (actualY == 0 && predY == 1)
        //  ++result[0][1];
        //else if (actualY == 1 && predY == 0)
        //  ++result[1][0];
        //else if (actualY == 1 && predY == 1)
        //  ++result[1][1];
      }
      return result;
    }

    public void ShowConfusionX(int[][] cm)
    {
      // simple display version
      Console.Write("actual = 0: ");
      Console.Write(cm[0][0].ToString().PadLeft(4) + " ");
      Console.WriteLine(cm[0][1].ToString().PadLeft(4));
      Console.Write("actual = 1: ");
      Console.Write(cm[1][0].ToString().PadLeft(4) + " ");
      Console.WriteLine(cm[1][1].ToString().PadLeft(4));
    }

    public void ShowConfusion(int[][] cm)
    {
      int n = cm.Length; // number classes
      double[] accuracies = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        int rowSum = 0;
        for (int j = 0; j "lt" n; ++j)
          rowSum += cm[i][j];
        accuracies[i] = (cm[i][i] * 1.0) / rowSum;
      }

      // assume binary classification
      int numCorrect = cm[0][0] + cm[1][1];
      int numItems = 
        cm[0][0] + cm[0][1] + cm[1][0] + cm[1][1];
      double overallAcc = (numCorrect * 1.0) / numItems;
     
      for (int i = 0; i "lt" n; ++i)
      {
        Console.Write(" actual " + 
          i.ToString().PadLeft(3) + " | ");
        for (int j = 0; j "lt" n; ++j)
          Console.Write(cm[i][j].ToString().PadLeft(6));
        Console.WriteLine("  | " + 
          accuracies[i].ToString("F4"));
      }
      Console.WriteLine("----------");
      Console.Write("predicted:    ");
      for (int j = 0; j "lt" cm[0].Length; ++j)
        Console.Write(j.ToString().PadLeft(6));
      Console.Write("  | " + overallAcc.ToString("F4"));
      Console.WriteLine("");
    }

  } // class LogisticRegression

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

} // ns

Training data:

# people_train_200.txt
# sex (0 = male, 1 = female)
# age, state (michigan, nebraska, oklahoma), income,
# politics type (conservative, moderate, liberal)
#
1, 0.24, 1, 0, 0, 0.2950, 0, 0, 1
0, 0.39, 0, 0, 1, 0.5120, 0, 1, 0
1, 0.63, 0, 1, 0, 0.7580, 1, 0, 0
0, 0.36, 1, 0, 0, 0.4450, 0, 1, 0
1, 0.27, 0, 1, 0, 0.2860, 0, 0, 1
1, 0.50, 0, 1, 0, 0.5650, 0, 1, 0
1, 0.50, 0, 0, 1, 0.5500, 0, 1, 0
0, 0.19, 0, 0, 1, 0.3270, 1, 0, 0
1, 0.22, 0, 1, 0, 0.2770, 0, 1, 0
0, 0.39, 0, 0, 1, 0.4710, 0, 0, 1
1, 0.34, 1, 0, 0, 0.3940, 0, 1, 0
0, 0.22, 1, 0, 0, 0.3350, 1, 0, 0
1, 0.35, 0, 0, 1, 0.3520, 0, 0, 1
0, 0.33, 0, 1, 0, 0.4640, 0, 1, 0
1, 0.45, 0, 1, 0, 0.5410, 0, 1, 0
1, 0.42, 0, 1, 0, 0.5070, 0, 1, 0
0, 0.33, 0, 1, 0, 0.4680, 0, 1, 0
1, 0.25, 0, 0, 1, 0.3000, 0, 1, 0
0, 0.31, 0, 1, 0, 0.4640, 1, 0, 0
1, 0.27, 1, 0, 0, 0.3250, 0, 0, 1
1, 0.48, 1, 0, 0, 0.5400, 0, 1, 0
0, 0.64, 0, 1, 0, 0.7130, 0, 0, 1
1, 0.61, 0, 1, 0, 0.7240, 1, 0, 0
1, 0.54, 0, 0, 1, 0.6100, 1, 0, 0
1, 0.29, 1, 0, 0, 0.3630, 1, 0, 0
1, 0.50, 0, 0, 1, 0.5500, 0, 1, 0
1, 0.55, 0, 0, 1, 0.6250, 1, 0, 0
1, 0.40, 1, 0, 0, 0.5240, 1, 0, 0
1, 0.22, 1, 0, 0, 0.2360, 0, 0, 1
1, 0.68, 0, 1, 0, 0.7840, 1, 0, 0
0, 0.60, 1, 0, 0, 0.7170, 0, 0, 1
0, 0.34, 0, 0, 1, 0.4650, 0, 1, 0
0, 0.25, 0, 0, 1, 0.3710, 1, 0, 0
0, 0.31, 0, 1, 0, 0.4890, 0, 1, 0
1, 0.43, 0, 0, 1, 0.4800, 0, 1, 0
1, 0.58, 0, 1, 0, 0.6540, 0, 0, 1
0, 0.55, 0, 1, 0, 0.6070, 0, 0, 1
0, 0.43, 0, 1, 0, 0.5110, 0, 1, 0
0, 0.43, 0, 0, 1, 0.5320, 0, 1, 0
0, 0.21, 1, 0, 0, 0.3720, 1, 0, 0
1, 0.55, 0, 0, 1, 0.6460, 1, 0, 0
1, 0.64, 0, 1, 0, 0.7480, 1, 0, 0
0, 0.41, 1, 0, 0, 0.5880, 0, 1, 0
1, 0.64, 0, 0, 1, 0.7270, 1, 0, 0
0, 0.56, 0, 0, 1, 0.6660, 0, 0, 1
1, 0.31, 0, 0, 1, 0.3600, 0, 1, 0
0, 0.65, 0, 0, 1, 0.7010, 0, 0, 1
1, 0.55, 0, 0, 1, 0.6430, 1, 0, 0
0, 0.25, 1, 0, 0, 0.4030, 1, 0, 0
1, 0.46, 0, 0, 1, 0.5100, 0, 1, 0
0, 0.36, 1, 0, 0, 0.5350, 1, 0, 0
1, 0.52, 0, 1, 0, 0.5810, 0, 1, 0
1, 0.61, 0, 0, 1, 0.6790, 1, 0, 0
1, 0.57, 0, 0, 1, 0.6570, 1, 0, 0
0, 0.46, 0, 1, 0, 0.5260, 0, 1, 0
0, 0.62, 1, 0, 0, 0.6680, 0, 0, 1
1, 0.55, 0, 0, 1, 0.6270, 1, 0, 0
0, 0.22, 0, 0, 1, 0.2770, 0, 1, 0
0, 0.50, 1, 0, 0, 0.6290, 1, 0, 0
0, 0.32, 0, 1, 0, 0.4180, 0, 1, 0
0, 0.21, 0, 0, 1, 0.3560, 1, 0, 0
1, 0.44, 0, 1, 0, 0.5200, 0, 1, 0
1, 0.46, 0, 1, 0, 0.5170, 0, 1, 0
1, 0.62, 0, 1, 0, 0.6970, 1, 0, 0
1, 0.57, 0, 1, 0, 0.6640, 1, 0, 0
0, 0.67, 0, 0, 1, 0.7580, 0, 0, 1
1, 0.29, 1, 0, 0, 0.3430, 0, 0, 1
1, 0.53, 1, 0, 0, 0.6010, 1, 0, 0
0, 0.44, 1, 0, 0, 0.5480, 0, 1, 0
1, 0.46, 0, 1, 0, 0.5230, 0, 1, 0
0, 0.20, 0, 1, 0, 0.3010, 0, 1, 0
0, 0.38, 1, 0, 0, 0.5350, 0, 1, 0
1, 0.50, 0, 1, 0, 0.5860, 0, 1, 0
1, 0.33, 0, 1, 0, 0.4250, 0, 1, 0
0, 0.33, 0, 1, 0, 0.3930, 0, 1, 0
1, 0.26, 0, 1, 0, 0.4040, 1, 0, 0
1, 0.58, 1, 0, 0, 0.7070, 1, 0, 0
1, 0.43, 0, 0, 1, 0.4800, 0, 1, 0
0, 0.46, 1, 0, 0, 0.6440, 1, 0, 0
1, 0.60, 1, 0, 0, 0.7170, 1, 0, 0
0, 0.42, 1, 0, 0, 0.4890, 0, 1, 0
0, 0.56, 0, 0, 1, 0.5640, 0, 0, 1
0, 0.62, 0, 1, 0, 0.6630, 0, 0, 1
0, 0.50, 1, 0, 0, 0.6480, 0, 1, 0
1, 0.47, 0, 0, 1, 0.5200, 0, 1, 0
0, 0.67, 0, 1, 0, 0.8040, 0, 0, 1
0, 0.40, 0, 0, 1, 0.5040, 0, 1, 0
1, 0.42, 0, 1, 0, 0.4840, 0, 1, 0
1, 0.64, 1, 0, 0, 0.7200, 1, 0, 0
0, 0.47, 1, 0, 0, 0.5870, 0, 0, 1
1, 0.45, 0, 1, 0, 0.5280, 0, 1, 0
0, 0.25, 0, 0, 1, 0.4090, 1, 0, 0
1, 0.38, 1, 0, 0, 0.4840, 1, 0, 0
1, 0.55, 0, 0, 1, 0.6000, 0, 1, 0
0, 0.44, 1, 0, 0, 0.6060, 0, 1, 0
1, 0.33, 1, 0, 0, 0.4100, 0, 1, 0
1, 0.34, 0, 0, 1, 0.3900, 0, 1, 0
1, 0.27, 0, 1, 0, 0.3370, 0, 0, 1
1, 0.32, 0, 1, 0, 0.4070, 0, 1, 0
1, 0.42, 0, 0, 1, 0.4700, 0, 1, 0
0, 0.24, 0, 0, 1, 0.4030, 1, 0, 0
1, 0.42, 0, 1, 0, 0.5030, 0, 1, 0
1, 0.25, 0, 0, 1, 0.2800, 0, 0, 1
1, 0.51, 0, 1, 0, 0.5800, 0, 1, 0
0, 0.55, 0, 1, 0, 0.6350, 0, 0, 1
1, 0.44, 1, 0, 0, 0.4780, 0, 0, 1
0, 0.18, 1, 0, 0, 0.3980, 1, 0, 0
0, 0.67, 0, 1, 0, 0.7160, 0, 0, 1
1, 0.45, 0, 0, 1, 0.5000, 0, 1, 0
1, 0.48, 1, 0, 0, 0.5580, 0, 1, 0
0, 0.25, 0, 1, 0, 0.3900, 0, 1, 0
0, 0.67, 1, 0, 0, 0.7830, 0, 1, 0
1, 0.37, 0, 0, 1, 0.4200, 0, 1, 0
0, 0.32, 1, 0, 0, 0.4270, 0, 1, 0
1, 0.48, 1, 0, 0, 0.5700, 0, 1, 0
0, 0.66, 0, 0, 1, 0.7500, 0, 0, 1
1, 0.61, 1, 0, 0, 0.7000, 1, 0, 0
0, 0.58, 0, 0, 1, 0.6890, 0, 1, 0
1, 0.19, 1, 0, 0, 0.2400, 0, 0, 1
1, 0.38, 0, 0, 1, 0.4300, 0, 1, 0
0, 0.27, 1, 0, 0, 0.3640, 0, 1, 0
1, 0.42, 1, 0, 0, 0.4800, 0, 1, 0
1, 0.60, 1, 0, 0, 0.7130, 1, 0, 0
0, 0.27, 0, 0, 1, 0.3480, 1, 0, 0
1, 0.29, 0, 1, 0, 0.3710, 1, 0, 0
0, 0.43, 1, 0, 0, 0.5670, 0, 1, 0
1, 0.48, 1, 0, 0, 0.5670, 0, 1, 0
1, 0.27, 0, 0, 1, 0.2940, 0, 0, 1
0, 0.44, 1, 0, 0, 0.5520, 1, 0, 0
1, 0.23, 0, 1, 0, 0.2630, 0, 0, 1
0, 0.36, 0, 1, 0, 0.5300, 0, 0, 1
1, 0.64, 0, 0, 1, 0.7250, 1, 0, 0
1, 0.29, 0, 0, 1, 0.3000, 0, 0, 1
0, 0.33, 1, 0, 0, 0.4930, 0, 1, 0
0, 0.66, 0, 1, 0, 0.7500, 0, 0, 1
0, 0.21, 0, 0, 1, 0.3430, 1, 0, 0
1, 0.27, 1, 0, 0, 0.3270, 0, 0, 1
1, 0.29, 1, 0, 0, 0.3180, 0, 0, 1
0, 0.31, 1, 0, 0, 0.4860, 0, 1, 0
1, 0.36, 0, 0, 1, 0.4100, 0, 1, 0
1, 0.49, 0, 1, 0, 0.5570, 0, 1, 0
0, 0.28, 1, 0, 0, 0.3840, 1, 0, 0
0, 0.43, 0, 0, 1, 0.5660, 0, 1, 0
0, 0.46, 0, 1, 0, 0.5880, 0, 1, 0
1, 0.57, 1, 0, 0, 0.6980, 1, 0, 0
0, 0.52, 0, 0, 1, 0.5940, 0, 1, 0
0, 0.31, 0, 0, 1, 0.4350, 0, 1, 0
0, 0.55, 1, 0, 0, 0.6200, 0, 0, 1
1, 0.50, 1, 0, 0, 0.5640, 0, 1, 0
1, 0.48, 0, 1, 0, 0.5590, 0, 1, 0
0, 0.22, 0, 0, 1, 0.3450, 1, 0, 0
1, 0.59, 0, 0, 1, 0.6670, 1, 0, 0
1, 0.34, 1, 0, 0, 0.4280, 0, 0, 1
0, 0.64, 1, 0, 0, 0.7720, 0, 0, 1
1, 0.29, 0, 0, 1, 0.3350, 0, 0, 1
0, 0.34, 0, 1, 0, 0.4320, 0, 1, 0
0, 0.61, 1, 0, 0, 0.7500, 0, 0, 1
1, 0.64, 0, 0, 1, 0.7110, 1, 0, 0
0, 0.29, 1, 0, 0, 0.4130, 1, 0, 0
1, 0.63, 0, 1, 0, 0.7060, 1, 0, 0
0, 0.29, 0, 1, 0, 0.4000, 1, 0, 0
0, 0.51, 1, 0, 0, 0.6270, 0, 1, 0
0, 0.24, 0, 0, 1, 0.3770, 1, 0, 0
1, 0.48, 0, 1, 0, 0.5750, 0, 1, 0
1, 0.18, 1, 0, 0, 0.2740, 1, 0, 0
1, 0.18, 1, 0, 0, 0.2030, 0, 0, 1
1, 0.33, 0, 1, 0, 0.3820, 0, 0, 1
0, 0.20, 0, 0, 1, 0.3480, 1, 0, 0
1, 0.29, 0, 0, 1, 0.3300, 0, 0, 1
0, 0.44, 0, 0, 1, 0.6300, 1, 0, 0
0, 0.65, 0, 0, 1, 0.8180, 1, 0, 0
0, 0.56, 1, 0, 0, 0.6370, 0, 0, 1
0, 0.52, 0, 0, 1, 0.5840, 0, 1, 0
0, 0.29, 0, 1, 0, 0.4860, 1, 0, 0
0, 0.47, 0, 1, 0, 0.5890, 0, 1, 0
1, 0.68, 1, 0, 0, 0.7260, 0, 0, 1
1, 0.31, 0, 0, 1, 0.3600, 0, 1, 0
1, 0.61, 0, 1, 0, 0.6250, 0, 0, 1
1, 0.19, 0, 1, 0, 0.2150, 0, 0, 1
1, 0.38, 0, 0, 1, 0.4300, 0, 1, 0
0, 0.26, 1, 0, 0, 0.4230, 1, 0, 0
1, 0.61, 0, 1, 0, 0.6740, 1, 0, 0
1, 0.40, 1, 0, 0, 0.4650, 0, 1, 0
0, 0.49, 1, 0, 0, 0.6520, 0, 1, 0
1, 0.56, 1, 0, 0, 0.6750, 1, 0, 0
0, 0.48, 0, 1, 0, 0.6600, 0, 1, 0
1, 0.52, 1, 0, 0, 0.5630, 0, 0, 1
0, 0.18, 1, 0, 0, 0.2980, 1, 0, 0
0, 0.56, 0, 0, 1, 0.5930, 0, 0, 1
0, 0.52, 0, 1, 0, 0.6440, 0, 1, 0
0, 0.18, 0, 1, 0, 0.2860, 0, 1, 0
0, 0.58, 1, 0, 0, 0.6620, 0, 0, 1
0, 0.39, 0, 1, 0, 0.5510, 0, 1, 0
0, 0.46, 1, 0, 0, 0.6290, 0, 1, 0
0, 0.40, 0, 1, 0, 0.4620, 0, 1, 0
0, 0.60, 1, 0, 0, 0.7270, 0, 0, 1
1, 0.36, 0, 1, 0, 0.4070, 0, 0, 1
1, 0.44, 1, 0, 0, 0.5230, 0, 1, 0
1, 0.28, 1, 0, 0, 0.3130, 0, 0, 1
1, 0.54, 0, 0, 1, 0.6260, 1, 0, 0

Test data:

# people_test_40.txt
#
0, 0.51, 1, 0, 0, 0.6120, 0, 1, 0
0, 0.32, 0, 1, 0, 0.4610, 0, 1, 0
1, 0.55, 1, 0, 0, 0.6270, 1, 0, 0
1, 0.25, 0, 0, 1, 0.2620, 0, 0, 1
1, 0.33, 0, 0, 1, 0.3730, 0, 0, 1
0, 0.29, 0, 1, 0, 0.4620, 1, 0, 0
1, 0.65, 1, 0, 0, 0.7270, 1, 0, 0
0, 0.43, 0, 1, 0, 0.5140, 0, 1, 0
0, 0.54, 0, 1, 0, 0.6480, 0, 0, 1
1, 0.61, 0, 1, 0, 0.7270, 1, 0, 0
1, 0.52, 0, 1, 0, 0.6360, 1, 0, 0
1, 0.30, 0, 1, 0, 0.3350, 0, 0, 1
1, 0.29, 1, 0, 0, 0.3140, 0, 0, 1
0, 0.47, 0, 0, 1, 0.5940, 0, 1, 0
1, 0.39, 0, 1, 0, 0.4780, 0, 1, 0
1, 0.47, 0, 0, 1, 0.5200, 0, 1, 0
0, 0.49, 1, 0, 0, 0.5860, 0, 1, 0
0, 0.63, 0, 0, 1, 0.6740, 0, 0, 1
0, 0.30, 1, 0, 0, 0.3920, 1, 0, 0
0, 0.61, 0, 0, 1, 0.6960, 0, 0, 1
0, 0.47, 0, 0, 1, 0.5870, 0, 1, 0
1, 0.30, 0, 0, 1, 0.3450, 0, 0, 1
0, 0.51, 0, 0, 1, 0.5800, 0, 1, 0
0, 0.24, 1, 0, 0, 0.3880, 0, 1, 0
0, 0.49, 1, 0, 0, 0.6450, 0, 1, 0
1, 0.66, 0, 0, 1, 0.7450, 1, 0, 0
0, 0.65, 1, 0, 0, 0.7690, 1, 0, 0
0, 0.46, 0, 1, 0, 0.5800, 1, 0, 0
0, 0.45, 0, 0, 1, 0.5180, 0, 1, 0
0, 0.47, 1, 0, 0, 0.6360, 1, 0, 0
0, 0.29, 1, 0, 0, 0.4480, 1, 0, 0
0, 0.57, 0, 0, 1, 0.6930, 0, 0, 1
0, 0.20, 1, 0, 0, 0.2870, 0, 0, 1
0, 0.35, 1, 0, 0, 0.4340, 0, 1, 0
0, 0.61, 0, 0, 1, 0.6700, 0, 0, 1
0, 0.31, 0, 0, 1, 0.3730, 0, 1, 0
1, 0.18, 1, 0, 0, 0.2080, 0, 0, 1
1, 0.26, 0, 0, 1, 0.2920, 0, 0, 1
0, 0.28, 1, 0, 0, 0.3640, 0, 0, 1
0, 0.59, 0, 0, 1, 0.6940, 0, 0, 1

Bottom line: Theory says the sophisticated, complex encoding technique should work better than a super simple encoding technique. In one experiment, the simple encoding technique worked slightly better.

I took a set of synthetic source training data that looks like this:

M,29,Austin,expert,68.5,timid,54600.00
F,51,Juneau,novice,60.0,brave,25700.00
M,32,Boston,master,73.5,alert,74900.00
M,19,Irvine,novice,61.0,brave,30300.00
. . .

Each row represents a person. The fields are sex, age, city (10 values), skill (3 values), height (inches), psychology (3 values), income (dollars). The goal is to predict income from the other six variables, using gradient boosting regression with decision tree weak learners.

There are 200 training items and 40 test items.

I normalized the numeric predictors, age and height, and the target variable income by dividing by 100, 100, 100,000 respectively. I did this just for convenience. The categorical predictors had to be converted to numbers because I wasn’t using a CART style tree.

There are two main schemes to encode categorical data for use with tree-based systems. The first encoding scheme is super simple: use ordinal encoding regardless of categorical type. So for sex, M = 0, F = 1. For city, Austin = 0, Boston = 1, Casper = 2, Denver = 3, Eugene = 4, Fresno = 5, Gallup = 6, Helena = 7, Irvine = 8, Juneau = 9. For skill, novice = 0, expert = 1, master = 2. For psychology, timid = 0, alert = 1, brave = 2.

The super simple normalized and encoded data looks like:

0, 0.29, 0, 1, 0.685, 0, 0.54600
1, 0.51, 9, 0, 0.600, 2, 0.25700
0, 0.32, 1, 2, 0.735, 1, 0.74900
0, 0.19, 8, 0, 0.610, 2, 0.30300
. . . 

With the super simple encoding, you can encode and then perform a train-test split. Easy.

The second encoding scheme is quite complicated. The complex normalized and encoded data looks like:

0, 0.29, 0.383784, 1, 0.685, 1, 0, 0, 0.54600
1, 0.51, 0.580719, 0, 0.600, 0, 0, 1, 0.25700
0, 0.32, 0.363510, 2, 0.735, 0, 1, 0, 0.74900
0, 0.19, 0.515110, 0, 0.610, 0, 0, 1, 0.30300
. . .

The city variable is smoothed-target encoded with smoothing parameter m=10. The psychology variable is one-hot encoded. When using smoothed target encoding, there is a risk of leakage. So you must first split, then compute the smoothed encodings on the training data, then apply to both train and test data.

I ran both versions of the encoded data through a from-scratch C# implementation of gradient boosting regression. The key output of the super simple encoding scheme:

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.8000

MSE train = 0.0000
MSE test = 0.0017

R2 train = 0.9994
R2 test = 0.9484

The key output of the complicated encoding scheme:

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.9360

The super simple encoding scheme is inferior to the complex encoding scheme, but the simple scheme worked slightly (but not statistically significant) better. The reasons why super simple encoding works well with gradient boosting regression are very complicated and subtle — quite literally pages of ideas and details.

I have seen this phenomenon many times in machine learning: complicated doesn’t always beat simple.

I have advanced college degrees in math and computer science, but my first degree was in psychology from UC Irvine. One of my requirements was a class in sociology. I remember learning a few things about strip club science.

Strip clubs and their equivalents have been around for literally thousands of years, dating back to ancient Greece and Rome.

I like to write code every day. Writing code is a skill that must be practiced, and besides, I enjoy writing code. Over the course of my career, I saw countless colleagues move from an engineering role to an engineering-management role, and they all lost their coding ability with surprising speed.

I decided to combine quadratic regression (QR) with gradient boosting regression (GBR). The idea is simple: create a container class that holds a QR model and a GBR model. After the two component models are trained independently, the hybrid model prediction is just a weighted average of the QR prediction and the GBR prediction.

For my demo, 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 values on each line are predictors. The last value on each line is the target to predict. There are 200 training items and 40 test items. The data was generated by a neural network with random weights and biases. Therefore, in theory, QR should do better than GBR.

The output of a demo, with the QR model prediction weight = 0.70 and the GBR prediction weight = 0.30:

Begin C# hybrid quadratic regression + 
 gradient boost regression

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 hybrid quadratic regression + 
 gradient boost regression model

Setting nTrees for GBR model = 200
Setting maxDepth GBR model = 3
Setting lrnRate GBR model = 0.1000
Setting QR predict weight = 0.70
Setting GBR predict weight = 0.30
Done

Training hybrid model
Done

QR model base weights:
 -0.2630  0.0354 -0.0421  0.0341 -0.1124

Evaluating hybrid model
Accuracy train (within 0.10) = 0.9300
Accuracy test (within 0.10) = 0.9250

MSE train = 0.0001
MSE test = 0.0004

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

Predicted y = 0.4839

End demo

The choice of combining QR and GBR is arbitrary. Any regression techniques can be used (for example, linear regression, kernel ridge regression, support vector regression, neural network regression, random forest regression, gradient boost regression). And you can use more than two regression techniques.

In practice, I like to use regression techniques that are much different from each other (because they can handle different types of data). As a practical matter, the fewer hyperparameters to tune the better.

With QR, I used a version that needs no parameters (no learning rate, no L2 regression, etc.) For gradient boosting regression, I used absolute minimal hyperparameters: nTrees, maxDepth, and lrnRate. I used sensible defaults for minSample (2), minLeaf (1), numSplitCols (-1 = all).

The Predict() and Train() methods are simple, because all of the work is done by the component models:

public class Hybid_QR_GBR
{
. . .
  public double Predict(double[] x)
  {
    double predQR = this.qrModel.Predict(x);
    double predGBR = this.gbrModel.Predict(x);
    double result = 
      (predQR * this.weights[0]) + 
      (predGBR * this.weights[1]);
    return result;
  }

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

  public void Train(double[][] trainX, double[] trainY)
  {
    this.qrModel.Train(trainX, trainY);
    this.gbrModel.Train(trainX, trainY);
  }
. . .

public class QuadraticRegressor { . . . }
public class QRHouseholder { . . . } // helper for QR
public class GradientBoostRegressor { . . . }
public class DecisionTreeRegressor { . . . } // help GBR

Instead of wrapping the hybrid regressor in a container class, Hybid_QR_NN_Regressor, it’s perfectly possible to just use the QuadraticRegressor class and the GradientBoostRegressor class separately, and then compute a weighted average of their Predict() methods. But putting everything into a container class keeps things more organized and easier to work with.

The scikit-learn library provides a VotingRegressor module that combines two or more scikit regression models.

Most of my friends, and me, get obsessed with all kinds of strange things. For reasons that are unknown to me, I am mildly obsessed with findling movies that feature an animal or 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.

In the “Legend of Ochi” (2025), a farm girl on a fictional rural island discovers a wounded baby creature (that chitters a lot). She learns its language and works to return it to its family, all while avoiding a dangerous hunting party. Everything about this movie should make me hate it — cute creature, tough adolescent girl, evil hunters, and so on. But I really liked this movie a lot. I give it a solid B+ rating.

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 Hybrid_QR_GBR
{
  internal class HybridProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# hybrid quadratic " +
        "regression + gradient boost regression ");

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

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

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

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

            // 2. create and train model
      Console.WriteLine("\nCreating hybrid quadratic " +
        "regression + gradient boost regression model ");
      int nTrees = 200;
      int maxDepth = 3;
      double lrnRate = 0.10;
      int seed = 0;
      double[] weights = new double[] { 0.70, 0.30 };

     Console.WriteLine("Setting nTrees for GBR model = " +
        nTrees);
      Console.WriteLine("Setting maxDepth GBR model = " +
        maxDepth);
      Console.WriteLine("Setting lrnRate GBR model = " +
        lrnRate.ToString("F4"));

      Console.WriteLine("Setting QR predict weight = " +
        weights[0].ToString("F2"));
      Console.WriteLine("Setting GBR predict weight = " +
        weights[1].ToString("F2"));

      Hybid_QR_GBR_Regressor model = 
        new Hybid_QR_GBR_Regressor(nTrees, maxDepth,
         lrnRate, seed, weights);
      Console.WriteLine("Done ");

      Console.WriteLine("\nTraining hybrid model ");
      model.Train(trainX, trainY); // train components           Console.WriteLine("Done ");

      // 3. show model weights
      Console.WriteLine("\nQR model base weights: ");
      int dim = trainX[0].Length;
      for (int i = 0; i "lt" dim; ++i)
        Console.Write(model.qrModel.weights[i].
          ToString("F4").PadLeft(8));
      Console.WriteLine("");
      // Console.WriteLine("\nModel quadratic weights: ");
      // Console.WriteLine("\nModel interaction weights: ");
      
      // 4. evaluate model
      Console.WriteLine("\nEvaluating hybrid model ");
      double accTrain = model.Accuracy(trainX, trainY, 0.10);
      Console.WriteLine("Accuracy 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"));

      // 5. use model
      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"));

      // 6. TODO: implement model Save() and Load()

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

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

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

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

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

  } // class Program

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

  public class Hybid_QR_GBR_Regressor
  {
    // training parameters for QR and GBR
    // QR: no parameters needed
    // GBR: nTrees, maxDepth, lrnRate.
    // uses minSample = 2, minLeaf = 1

    public int nNeighbors;
    private Random rnd;  // not used w/ Pinv training
    public double[] weights = new double[0]; // each model

    public QuadraticRegressor qrModel;
    public GradientBoostRegressor gbrModel;

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

    public Hybid_QR_GBR_Regressor(int nTrees, int maxDepth,
      double lrnRate, int seed, double[] weights)
    {
      this.nNeighbors = nNeighbors;
      this.rnd = new Random(seed);
      this.weights = weights;

      this.qrModel =
        new QuadraticRegressor(seed);
      this.gbrModel =
        new GradientBoostRegressor(nTrees: nTrees, 
        maxDepth: maxDepth, minSamples: 2,
        minLeaf: 1, numSplitCols: -1, lrnRate: lrnRate,
        seed);
    }

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

    public double Predict(double[] x)
    {
      double predQR = this.qrModel.Predict(x);
      double predGBR = this.gbrModel.Predict(x);
      double result = (predQR * this.weights[0]) +
        (predGBR * this.weights[1]);
      return result;
    }

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

    public void Train(double[][] trainX, double[] trainY)
    {
      this.qrModel.Train(trainX, trainY);
      this.gbrModel.Train(trainX, trainY);
    }

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

    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 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) "lte"
          Math.Abs(pctClose * actualY))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

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

  } // class Hybid_QR_NN


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

  public class QuadraticRegressor
  {
    public double[] weights;  // regular, quad, interactions
    public double bias;
    private Random rnd;  // not used w/ Pinv training

    public QuadraticRegressor(int seed = 0)
    {
      this.weights = new double[0];  // empty, but not null
      this.bias = 0; // dummy value
      this.rnd = new Random(seed);
    }

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

    public double Predict(double[] x)
    {
      int dim = x.Length;
      double result = 0.0;

      int p = 0; // points into this.weights
      for (int i = 0; i "lt" dim; ++i)   // base terms
        result += x[i] * this.weights[p++];

      for (int i = 0; i "lt" dim; ++i)  // quadratic terms
        result += (x[i] * x[i]) * this.weights[p++];

      for (int i = 0; i "lt" dim - 1; ++i)  // interactions
        for (int j = i + 1; j "lt" dim; ++j)
          result += (x[i] * x[j]) * this.weights[p++];

      result += this.bias;
      return result;
    }

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

    public void Train(double[][] trainX, double[] trainY)
    {
      // train using MP pseudo-inverse QR-Householder
      // no regulaization
      // w = pinv(designX) * y
      int nRows = trainX.Length; // not used
      int dim = trainX[0].Length;
      int nInteractions = (dim * (dim - 1)) / 2;
      this.weights = new double[dim + dim + nInteractions];

      double[][] Xa = MatAugment(trainX);  // add quad cols
      double[][] X = MatToDesign(Xa);  // add 1.0s col

      double[][] Xpinv = QRHouseholder.MatPinv(X);

      double[] biasAndWts = MatVecProd(Xpinv, trainY);
      this.bias = biasAndWts[0];  // bias is at [0]
      for (int i = 1; i "lt" biasAndWts.Length; ++i)
        this.weights[i - 1] = biasAndWts[i];
      return;
    }

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

    private static double[][] MatAugment(double[][] trainX)
    {
      // add quadratic and interaction columns
      int nRows = trainX.Length;  // src and dest
      int dim = trainX[0].Length;  // src
      int nInteractions = dim * (dim - 1) / 2;
      int nColsDest = dim + dim + nInteractions;

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

      for (int i = 0; i "lt" nRows; ++i)
      {
        int p = 0; // points to column of result

        for (int j = 0; j "lt" dim; ++j) // base
          result[i][p++] = trainX[i][j];

        for (int j = 0; j "lt" dim; ++j) // quadratic
          result[i][p++] = trainX[i][j] * trainX[i][j];

        for (int j = 0; j "lt" nInteractions - 1; ++j)
          for (int k = j + 1; k "lt" dim; ++k)
            result[i][p++] = trainX[i][j] * trainX[i][k];
      }

      return result;
    }

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

    private static double[][] MatToDesign(double[][] X)
    {
      // add leading column of 1.0s to handle bias term
      int nRows = X.Length;  // src and dest
      int dim = X[0].Length;

      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[dim + 1]; // extra col

      for (int i = 0; i "lt" nRows; ++i)
      {
        result[i][0] = 1.0;
        for (int j = 1; j "lt" result[0].Length; ++j)
        {
          result[i][j] = X[i][j - 1];
        }
      }
      return result;
    }

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

    private static double[] MatVecProd(double[][] M,
      double[] v)
    {
      // helper for Train()
      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;
    }

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

  } // class QuadraticRegressor

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

  public class QRHouseholder
  {
    // for class QuadraticRegressor Train()
    // container for MP pseudo-inverse via QR-Householder
    // A = Q * R
    // pinv(A) = inv(R) * inv(Q)  note order matters
    //         = inv upper tri (easy) * transpose (easy)

    public static double[][] MatPinv(double[][] M)
    {
      double[][] Q; double[][] R;
      MatDecompQR(M, out Q, out R);  // Householder
      double[][] Ri = MatInvUpperTri(R);
      double[][] Qi = MatTranspose(Q);
      double[][] result = MatProduct(Ri, Qi);
      return result;
    }

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

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

      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      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];
        }
      }
      return result;
    }

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

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

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

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

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

    public 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 = new double[aRows][];
      for (int i = 0; i "lt" aRows; ++i)
        result[i] = new double[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;
    }

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

    public static void MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // Householder algorithm
      int m = A.Length; int n = A[0].Length;
      if (m "lt" n)
        Console.WriteLine("FATAL: nRows must be gte nCols");

      double[][] QQ = MatMake(m, m); // working full Q
      for (int i = 0; i "lt" m; ++i)
        QQ[i][i] = 1.0;  // identity matrix

      double[][] RR = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          RR[i][j] = A[i][j]; // copy of A is working R

      int k = Math.Min(m, n);  // or just use n
      for (int j = 0; j "lt" k; ++j) // main processing loop
      {
        int xn = m - j;
        double[] x = new double[xn];
        for (int i = 0; i "lt" xn; ++i)
          x[i] = RR[j + i][j];

        double ss = 0.0;
        for (int i = 0; i "lt" xn; ++i)
          ss += x[i] * x[i];
        double normX = Math.Sqrt(ss);

        // if (normX == 0.0) continue;  // risky
        if (Math.Abs(normX) "lt" 1.0e-12) continue;

        double sign;
        if (x[0] "gte" 0.0) sign = -1.0;
        else sign = 1.0; // counter-intuitive

        double[] u = new double[xn];
        for (int i = 0; i "lt" xn; ++i)
          u[i] = x[i] / (x[0] - sign * normX); // check div 0
        u[0] = 1.0;

        // compute scaling factor tau = 2 / (u^T * u)
        double tau = -sign * (x[0] - sign * normX) / normX;

        // dimensions for sub-matrices
        int nRowsSubR = m - j; int nColsSubR = n - j;
        int nRowsSubQ = m; int nColsSubQ = m - j;

        double[] vr = new double[nColsSubR];
        for (int c = 0; c "lt" nColsSubR; ++c)
        {
          double acc = 0.0;
          for (int r = 0; r "lt" nRowsSubR; ++r)
            acc += u[r] * RR[j + r][j + c];
          vr[c] = acc;
        }

        double[] vq = new double[nRowsSubQ];
        for (int r = 0; r "lt" nRowsSubQ; ++r)
        {
          double acc = 0.0;
          for (int c = 0; c "lt" nColsSubQ; ++c)
            acc += u[c] * QQ[r][j + c];
          vq[r] = acc;
        }

        // update sub-R
        for (int r = 0; r "lt" nRowsSubR; ++r)
          for (int c = 0; c "lt" nColsSubR; ++c)
            RR[j + r][j + c] -= tau * u[r] * vr[c];

        // update sub-Q
        for (int r = 0; r "lt" nRowsSubQ; ++r)
          for (int c = 0; c "lt" nColsSubQ; ++c)
            QQ[r][j + c] -= tau * vq[r] * u[c];

      } // j

      // extract QQ RR into out params
      Q = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          Q[i][j] = QQ[i][j];

      R = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          R[i][j] = RR[i][j];

      return;
    } // MatDecompQR

  } // class QRHouseholder

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

  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
      {
        // compute residuals from prev tree predictions
        for (int i = 0; i "lt" n; ++i)
          residuals[i] = trainY[i] - preds[i];

        int treeSeed = this.rnd.Next();
        DecisionTreeRegressor dtr =
          new DecisionTreeRegressor(
            this.maxDepth,
            this.minSamples,
            this.minLeaf,
            this.numSplitCols,
            false, // don't save rows info
            treeSeed);
        dtr.Train(trainX, residuals); // predict residuals

        // predict the 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)
    {
      double result = this.pred0;
      for (int t = 0; t "lt" this.nTrees; ++t)
      {
        double predResidual = this.trees[t].Predict(x);
        double delta = this.lrnRate * predResidual;
        result += delta;
      }
      return result;
    }

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

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

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

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

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

I won’t try to explain SVR or SMO. Both topics are extremely complicated. But I will caution you that I’ve seen a huge amount of grossly incorrect information about these two topics on the Internet.

Here’s the output of my demo:

Begin C# support vector regression with SMO

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

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

First three target y:
  0.4840
  0.1568
  0.8054

Creating SVR object
Setting RBF gamma = 0.30
Setting epsilon = 0.0320
Setting C = 1.00
Setting KKT tol = 0.000010
Setting maxPasses = 100
Done

Training SVR model using SMO
Done

Model dual weights:
 -0.0832  -0.2560   0.4079  -0.1501 . . .   0.0578
  0.0903  -0.0010  -0.0065   0.2746 . . .   0.1430
. . .
  0.0855  -0.2647   0.0144   0.3009 . . .  -0.1137
  0.2584

Model bias = 0.4156

Number supp vectors = 196

Evaluating model

Train acc (within 0.10) = 0.9200
Test acc (within 0.10) = 0.9250

Train MSE = 0.0002
Test MSE = 0.0002

Predicting for trainX[0]
Predicted y = 0.4811

End SVR with SMO training demo

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

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

Support vector regression with quadratic programming training (not SMO). But illustrates epsilon and C. The slack variables are used in QP training.

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

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

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

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

Using SVR with SMO training involves specifying 5 hyperparameters: the internal RBF gamma, the epsilon “closeness” value, C regularization, KKT tolerance, and maxPasses. SVR is extremely sensitive to these hyperparameter values, and they interact in complex ways, which makes training and tuning an SVR-SMO model extremely difficult.

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

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.

So I better keep my day job as the saying goes.

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

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

// SVR using sequential minimal optimization training (SMO)
// SVR is extrememly prone to model overfitting
// (libsvm library is very complicated)

namespace SupportVectorRegressionSMO
{
  internal class SupportVectorRegressionSMOProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# support " +
        "vector regression with SMO ");

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

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

      Console.WriteLine("\nFirst three X predictors: ");
      for (int i = 0; i "lt" 3; ++i)
        Utils.VecShow(trainX[i], 4, 9);
      Console.WriteLine("\nFirst three target y: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine(trainY[i].ToString("F4").
          PadLeft(8));

      Console.WriteLine("\nCreating SVR object");
      double gamma = 0.30;    // RBF param
      double epsilon = 0.032;
      double C = 1.0;
      double tol = 1.0e-5;  // KKT tolerance
      int maxPasses = 100;  // number times no change
      
      Console.WriteLine("Setting RBF gamma = " +
        gamma.ToString("F2"));
      Console.WriteLine("Setting epsilon = " +
        epsilon.ToString("F4"));
      Console.WriteLine("Setting C = " +
        C.ToString("F2"));
      Console.WriteLine("Setting KKT tol = " +
        tol.ToString("F6"));
      Console.WriteLine("Setting maxPasses = " +
        maxPasses);
      
      SVR model = new SVR(gamma, epsilon, C, tol,
        maxPasses, seed: 0);
      Console.WriteLine("Done ");

      Console.WriteLine("\nTraining SVR model using" +
        " SMO ");
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");
      
      //Console.WriteLine("\nModel alpha: ");
      //Utils.VecShow(model.alpha, 4, 9);
      //Console.WriteLine("\nModel alpha*: ");
      //Utils.VecShow(model.alphaStar, 4, 9);

      Console.WriteLine("\nModel dual weights: ");
      Utils.VecShow(model.dualWeights, 4, 9);
      Console.WriteLine("\nModel bias = " + 
        model.b.ToString("F4"));
      Console.WriteLine("\nNumber supp vectors = " + 
        model.dualWeights.Length);

      Console.WriteLine("\nEvaluating model ");
      double trainAcc =
        model.Accuracy(trainX, trainY, 0.10);
      double testAcc =
        model.Accuracy(testX, testY, 0.10);

      Console.WriteLine("\nTrain acc (within 0.10) = " +
        trainAcc.ToString("F4"));
      Console.WriteLine("Test acc (within 0.10) = " +
        testAcc.ToString("F4"));

      double trainMSE = model.MSE(trainX, trainY);
      double testMSE = model.MSE(testX, testY);

      Console.WriteLine("\nTrain MSE = " +
        trainMSE.ToString("F4"));
      Console.WriteLine("Test MSE = " +
        testMSE.ToString("F4"));

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

      Console.WriteLine("\nEnd SVR with SMO training demo ");
      Console.ReadLine();
    } // Main()

  } // class Program

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

  public class SVR
  {
    public double gamma;  // for RBF kernel
    public double epsilon;
    public double C; // regularization
    public double tol;  // for KKT conditions
    public int maxPasses;  // times no improvement
    public double[][] suppX;
    public double[] suppY;
    public double[] alpha;
    public double[] alphaStar;
    public double[] dualWeights; // alpha* - alpha
    public double b; // bias
    public Random rnd;  // training is probabilistic

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

    public SVR(double gamma, double epsilon, double C, 
      double tol, int maxPasses, int seed = 0)
    {
      this.gamma = gamma;
      this.epsilon = epsilon;
      this.C = C;
      this.tol = tol;
      this.maxPasses = maxPasses;

      this.suppX = new double[0][]; // compiler happy
      this.suppY = new double[0];
      this.alpha = new double[0];
      this.alphaStar = new double[0];
      this.dualWeights = new double[0];
      this.b = 0.0;
      this.rnd = new Random(seed); // SMO selection
    } // ctor

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

    public void Train(double[][] trainX, double[] trainY)
    {
      // train using SMO algorithm
      this.suppX = trainX;  // needed by Predict()
      this.suppY = trainY;  // not used this version
      int n = trainX.Length;

      this.alpha = new double[n];
      this.alphaStar = new double[n];
      this.b = Utils.VecMean(trainY);

      // precompute all rbf values to K for fast train
      double[][] K = this.MakeK(trainX);

      int nPasses = 0;  // number times no change in alpha
      while (nPasses "lt" this.maxPasses)
      {
        int numChangedAlphas = 0;

        for (int i = 0; i "lt" n; ++i) // walk each train
        {
          double[] iCurrCol = Utils.MatGetColumn(K, i);
          double iPred = 
            Utils.VecDot(Utils.VecSubtract(this.alphaStar,
            this.alpha), iCurrCol) + this.b;
          double iErr = iPred - trainY[i];

          // check KKT to see if [i] warrants update
          if ((iErr "gt" this.epsilon + this.tol &&
            this.alpha[i] "lt" this.C) ||
            (iErr "gt" this.epsilon + this.tol &&
            this.alphaStar[i] "gt" 0.0) ||
            (iErr "lt" -this.epsilon - this.tol &&
            this.alpha[i] "gt" 0.0) ||
            (iErr "lt" -this.epsilon - this.tol &&
            this.alphaStar[i] "lt" this.C))
          {
            // update to [i] is warranted
            // pick random j, not the same as i
            int j = i;
            while (j == i)
              j = this.rnd.Next(0, n);

            double[] jCurrCol = Utils.MatGetColumn(K, j);
            double jPred = 
              Utils.VecDot(Utils.VecSubtract(this.alphaStar,
              this.alpha), jCurrCol) + this.b;
            double jErr = jPred - trainY[j];

            // save old alpha and alphaStar
            double iOldAlpha = this.alpha[i];
            double iOldAlphaStar = this.alphaStar[i];
            double jOldAlpha = this.alpha[j];
            double jOldAlphaStar = this.alphaStar[j];

            // kernel second step denominator
            double eta = K[i][i] + K[j][j] - (2.0 * K[i][j]);
            if (eta "lte" 0.0)
              continue;

            // SVR linear constraint constant
            double constraint = 
              (iOldAlphaStar - iOldAlpha) +
              (jOldAlphaStar - jOldAlpha);

            // joint proxy variable
            double jOldDelta = 
              jOldAlphaStar - jOldAlpha;  // jOldDelta
            double jNewDelta = 
              jOldDelta + (iErr - jErr) / eta; // jNewDelta

            // compute bounds
            double L = Math.Max(-this.C,
              constraint - this.C);
            double H = Math.Min(this.C,
              constraint + this.C);
            // clip
            if (jNewDelta "lt" L)
              jNewDelta = L;
            else if (jNewDelta "gt" H)
              jNewDelta = H;

            if (Math.Abs(jNewDelta - jOldDelta) "lt"
              1.0e-5)
                continue;
            
            // reconstruct j
            if (jNewDelta "gte" 0.0)
            {
              this.alphaStar[j] = jNewDelta;
              this.alpha[j] = 0.0;
            }
            else
            {
              this.alphaStar[j] = 0.0;
              this.alpha[j] = -jNewDelta;
            }

            // update i
            double iNewValue = 
              constraint - jNewDelta;  // iNewValue
            if (iNewValue "gte" 0.0)
            {
              this.alphaStar[i] = iNewValue;
              this.alpha[i] = 0.0;
            }
            else
            {
              this.alphaStar[i] = 0.0;
              this.alpha[i] = -iNewValue;
            }

            // update bias
            double b1 = this.b - iErr - ((this.alphaStar[i] -
              this.alpha[i]) - (iOldAlphaStar - iOldAlpha)) *
              K[i][i] - ((this.alphaStar[j] -
              this.alpha[j]) - (jOldAlphaStar - jOldAlpha)) *
              K[i][j];


            double b2 = this.b - jErr - ((this.alphaStar[i] -
              this.alpha[i]) - (iOldAlphaStar - iOldAlpha)) *
              K[i][j] - ((this.alphaStar[j] -
              this.alpha[j]) - (jOldAlphaStar - jOldAlpha)) *
              K[i][j];

            if (0.0 "lt" this.alpha[i] &&
                this.alpha[i] "lt" this.C ||
                0.0 "lt" this.alphaStar[i] && 
                this.alphaStar[i] "lt" this.C)
              this.b = b1;
            else if (0.0 "lt" this.alpha[j] &&
                this.alpha[j] "lt" this.C ||
                0.0 "lt" this.alphaStar[j] &&
                this.alphaStar[j] "lt" this.C)
              this.b = b2;
            else
              this.b = (b1 + b2) / 2.0;

            ++numChangedAlphas;
          }
          else
          {
            ; // no change to alpha[i] warranted
          }
          
        } // each item i

        if (numChangedAlphas == 0)
          ++nPasses;
        else
          nPasses = 0;  // reset
      } // while

      // prune
      // 1. combine alpha and alpha*
      this.dualWeights = Utils.VecSubtract(this.alphaStar,
        this.alpha);

      // 2. compute mask based on small dual wts
      List"lt"int"gt" lstMask = new List"lt"int"gt"();
      double maskTol = 1.0e-4;  // larger tol == fewer svs
      for (int i = 0; i "lt" n; ++i)
        if (Math.Abs(this.dualWeights[i]) "gt" maskTol)
          lstMask.Add(i);
      int[] svMask = lstMask.ToArray();

      // 3. mask supp vecs
      this.suppX = Utils.MatSelectRows(trainX, svMask);
      this.suppY = Utils.VecSelectItems(trainY, svMask);

      // 4. mask dual weights
      this.dualWeights = 
        Utils.VecSelectItems(this.dualWeights, svMask);

      // leave source alphas alone . . 

      return;  // all done
    } // Train

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

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

    private double[][] MakeK(double[][] X)
    {
      // Kernel-Gram matrix helper for Train()
      // pre-compute all similarities, to avoid re-computes
      int n = X.Length;
      double[][] result = Utils.MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = this.RBF(this.suppX[i],
          this.suppX[j]);
      return result;
    }

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

    public double Predict(double[] x)
    {
      int n = this.suppX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
      {
        double k = this.RBF(this.suppX[i], x);
        //sum += (this.alphaStar[i] - this.alpha[i]) * k;
        sum += this.dualWeights[i] * k;
      }
      return sum + this.b;
    }

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

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

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

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

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

  public class Utils
  {
    // ------------------------------------------------------

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

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

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

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

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

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

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

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

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

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

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

    public static double VecMean(double[] vec)
    {
      int n = vec.Length;  // not 0
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += vec[i];
      return sum / n; 
    }

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

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

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

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

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

    public static double[] VecSelectItems(double[] vec,
      int[] idxs)
    {
      int n = idxs.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        result[i] = vec[idxs[i]];
      }
      return result;
    }

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

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

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

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

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

    public static 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 Utils

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

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

A simple decision tree for regression will almost always overfit the training data where the model accuracy on the training data is near 100% but accuracy on new, previously unseen data is very poor.

The idea of AdaBoost.R2 (adaptive boosting for regression, version 2) regression is to create a collection/ensemble of many simple decision trees, where each tree is trained on a different subset of the source training data. Each subset of training data is selected to emphasize items that were incorrectly predicted by the previous tree.

As each tree is trained, it receives an importance weight (sometimes labeled as “alpha”) where a higher importance weight means a better predictor. Then, after all the trees are trained, a prediction for an input vector x is the weighted (by tree importance) median of the predictions of the trees in the collection.

A much simpler tree-based regression technique, called bagging (bootstrap aggregation) tree regression, was published in 1994. Bagging tree regression evolved into random forest regression (1995), which evolved into AdaBoost.R and AdaBoost.R2 regression (1997), which evolved into gradient boosting regression (1999) — a fascinating (but very long) story. Maybe I’ll post that story some day.

I used AI to optimize my AdaBoost.R2 regression system. AI found many optimizations so I added those optimizations and I made a demo. The demo data is synthetic. It was generated by a 5-10-1 neural network with random weight and bias values. 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
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
. . .

The first five values on each line are predictors. The last value on each line is the target y value to predict. There are 200 training items and 40 test items. Here’s the output of a demo:

Begin AdaBoost.R2 (tree) regression from scratch C#

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

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

First three train y:
  0.4840
  0.1568
  0.8054

Setting nEstimators = 300
Setting lrnRate = 1.0000
Setting tree maxDepth = 4
Setting tree minSamples = 2
Setting tree minLeaf = 1
Setting tree numSplitCols = -1

Training AdaBoost.R2 model
Done
Created 300 estimators

Accuracy train (within 0.10): 0.7100
Accuracy test (within 0.10): 0.5500

MSE train: 0.0006
MSE test: 0.0014

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

End demo

Notice that the model still overfits significantly. This is because the synthetic data is more or less spread evenly, which is not a good use-case for AdaBoost.

Because of this weakness, the more powerful gradient boosting regression model is frequently used instead of AdaBoost.R2 regression or bagging tree regression or random forest regression. That said however, for some problem scenarios, AdaBoost.R2 regression work quite well.

Machine learning regression algorithms don’t have any special super powers — they’re just clever mathematics. When I was young I learned to read from comic books. One of my favorite heroes had no special powers. Green Arrow, and his sidekick Speedy, used their wits, amazing archery skill, and clever trick arrows to defeat villains.

Right: “Brave and Bold” #50 (October-November 1963) was an anthology title that paired up a couple of heroes. This issue featured Green Arrow with the Martian Manhunter. Cover art by George Roussos.

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

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

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

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

      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        new int[] { colY }, ',', "#"));

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

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

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

      // 2. create and train model
      int nEstimators = 300;
      double lrnRate = 0.90; // indirect regularizer
      int maxDepth = 4;
      int minSamples = 2; // using minLeaf = 1

      // The original AdaBoost.R2 algorithm paper does
      // not use a lrn_rate. The original item weight
      // update is:
      //   w = w * beta^(1-L), beta = L / (1-L)
      // where L is the loss and beta is a confidence.
      // The scikit AdaBoostRegressor adds a lrn_rate:
      //   w = w * (beta^(1-L))^lrn_rate, beta = L / (1-L)
      // If lrn_rate = 1.0, the modified form is the
      // same as the original form. Smaller lrn_rate
      // less than 1.0 gives smaller updates.

      // In AdaBoost.R2, because predictions are combined
      // via a weighted median of trees predicting full
      // y-values, tweaking the learning rate alters how
      // aggressively the algorithm shifts focus to
      // difficult data points, rather than directly
      // shrinking prediction steps. This makes tuning
      // learning rate in AdaBoost feel much less like
      // smooth "step-size scaling" and more like adjusting
      // the model's memory/aggressiveness.

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

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

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

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

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

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

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

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

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

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

  } // class Program

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

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

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

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

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

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

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

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

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

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

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

        // 3. train base tree on the bootstrap subset
        int treeSeed = this.rnd.Next(0, 1_000_000);
        DecisionTreeRegressor dtr = 
          new DecisionTreeRegressor(
          maxDepth: this.maxDepth,
          minSamples: this.minSamples,
          minLeaf: 1,
          numSplitCols: -1,
          saveRows: false, // don't save rows info
          seed: treeSeed
        );
        dtr.Train(subsetX, ySubset);

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

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

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

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

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

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

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

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

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

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

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

      int[] result = new int[size];

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

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

      return result;
    } // MyChoice()

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

    private static int SearchCdf(double[] cdf, double target)
    {
      // binary search to isolate the target interval
      int low = 0;
      int high = cdf.Length - 1;
            
      while (low "lte" high)
      {
        int mid = low + (high - low) / 2;
        if (cdf[mid] "gte" target)
          high = mid - 1;
        else
          low = mid + 1;
      }

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

      return low;
    }

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

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

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

    public double Predict(double[] x)
    {
      int nTrees = this.estimators.Count;
      double[] preds = new double[nTrees];
      double[] modelWts = new double[nTrees];

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

      return WeightedMedian(preds, modelWts);
    }

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

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

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

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

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

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

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

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

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

  } // class AdaBoostRegressor

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

  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

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


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

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