James D. McCaffrey

There are three ways to train a basic linear regression model: 1.) using stochastic gradient descent, 2.) using left pseudo-inverse (normal equations) via Cholesky inverse, 3.) using relaxed Moore-Penrose pseudo-inverse via one of many possible inverses.

I have implemented many of the training algorithms from scratch. For type (3) training, my approach-of-choice is QR-Householder. One of the (very) minor disadvantages of using MP pseudo-inverse training is that there’s no easy way to add L2 regularization. Before I go any further, let me state forcefully that I almost never use regularization with linear regression. It’s just not useful, but that’s another (surprisingly complicated) story.

One evening, after I walked my dogs, just for fun, I figured I’d add L2 regularization to my implementation of linear regression trained with MP pseudo-inverse via QR-Householder.

Output of the demo:

Begin C# linear regression using MP pinv
 (QR-Householder) training with L2 regularization

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

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

First three train y:
  0.4840
  0.1568
  0.8054

Creating and training Linear Regression model
 using QR p-inverse
Setting L2 lamda = 1.0000
Done

Coefficients/weights:
-0.2618  0.0331  -0.0453  0.0353  -0.1132
Bias/constant: 0.3618

Evaluating model

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

MSE train = 0.0026
MSE test = 0.0019

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

Predicted y = 0.5311

End demo

The demo data is synthetic. It was generated by a neural network, and so linear regression cannot predict it very well. There are 200 training items and 40 test items so it’s a small dataset. There are 5 predictors.

Adding L2 regularization to a linear regression model that’s trained using MP pseudo-inverse is tricky. Expressed in a diagram:

The diagram illustrates an example with just 4 predictor variables and just 20 training items. First the training data is converted to a design matrix by adding a leading column of 1.0 values. Next a special regularization matrix is added to the bottom of the design matrix. The square root of the regularization constant is added to the diagonal of the regularization matrix, except for the upper-left cell which corresponds to the model bias. The target y values have associated 0.0s appended as shown.

At this point the modified X matrix and y vector are trained using MP pinv as usual.

If you get the feeling that this is complicated, well it is, especially since regularization for linear regression is rarely useful.

Interesting experiment but not really useful in a practical way. The technique can be used with any linear model, notably quadratic regression where it can sometimes be mildly useful

I am fascinated to the point of obsession by mathematics, computer science, and machine learning. All of these things represent a model of reality in some way.

When I was a young man, I loved building plastic models, especially those from the Revell company. Some of my all time favorites were models produced in the 1950s that were space related. Here’s one model that I remember with great fondness.

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

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

namespace LinearRegressionPinvQRHouseholder
{
  internal class LinearRegressionPinvProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression" +
        " using MP pinv (QR-Householder) training with" +
        " L2 regularization ");

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

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

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

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

      // 2. create and train model using pseudo-inverse
      Console.WriteLine("\nCreating and training" +
        " Linear Regression model using QR p-inverse");
      LinearRegressor model = new LinearRegressor();

      // train no regularization
      //model.Train(trainX, trainY);
      //Console.WriteLine("Done ");

      // train with regularization
      double lamda = 1.0;
      Console.WriteLine("Setting L2 lamda = " +
        lamda.ToString("F4"));
      model.Train(trainX, trainY, lamda);
      Console.WriteLine("Done ");

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

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

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

      // 4. use model to predict first training item
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 9);
      double predY = model.Predict(x);
      Console.WriteLine("\nPredicted y = " +
        predY.ToString("F4"));

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

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

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

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

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

  } // class Program

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

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

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

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

    // ------------------------------------------------------
    // primary: Train(), Predict(), Accuracy(), MSE()
    // helpers: MatToDesign(), MatVecProd()
    // ------------------------------------------------------

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

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

    public void Train(double[][] trainX, double[] trainY)
    {
      // no regularization.
      // wts = pinv(designX) * y
      int dim = trainX[0].Length;
      this.weights = new double[dim];
      double[][] X = MatToDesign(trainX);
      double[][] Xpinv = QRHouseholder.MatPinv(X);
      double[] biasAndWts = MatVecProd(Xpinv, trainY);
      this.bias = biasAndWts[0];
      for (int i = 1; i "lt" biasAndWts.Length; ++i)
        this.weights[i - 1] = biasAndWts[i];
      return;
    }

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

    public void Train(double[][] trainX, double[] trainY,
      double lamda)
    {
      // train using MP pinv QR Householder with L2
      int nRows = trainX.Length;  // 200
      int dim = trainX[0].Length;  // 5
      this.weights = new double[dim];

      double[][] Xd = MatToDesign(trainX); // 200 by 6
      double[][] X = MatRegularize(Xd, lamda); // tricky
      
      double[] Y = new double[nRows + dim + 1];
      for (int j = 0; j "lt" nRows; ++j)
        Y[j] = trainY[j];

      double[][] Xpinv = QRHouseholder.MatPinv(X);
      double[] biasAndWts = MatVecProd(Xpinv, Y);
      this.bias = biasAndWts[0];
      for (int i = 1; i "lt" biasAndWts.Length; ++i)
        this.weights[i - 1] = biasAndWts[i];
      return;
    }

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

    private static double[][] MatRegularize(double[][] Xd,
      double lamda)
    {
      // Xd is a design matrix with leading col of 1.0s
      // add dim by dim to bottom of X
      // then . . . 
      int nRows = Xd.Length;  // src
      int nCols = Xd[0].Length;  // src
      double[][] result = MatMake(nRows + nCols, nCols);

      // copy top into result, inc leading 1.0s
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = Xd[i][j];
        
      // fill bottom starting at [1][1] (skip bias)
      int col = 1;
      for (int i = nRows + 1; i "lt" result.Length; ++i)
        result[i][col++] = Math.Sqrt(lamda);
      
      return result;
    }

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

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

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

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

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

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

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

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

      return result;
    }

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

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

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

  } // class LinearRegressor

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

  public class QRHouseholder
  {
    // 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)
    {
      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;
        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


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

} // ns

Training data:

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

Test data:

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

I came across an anomaly detection idea that I hadn’t seen before. It’s called an Isolation Forest. Briefly, if you send data to a decision tree regressor, data items that are anomalous will end up closer to the root node of the tree.

For example, suppose you have just four data items:

[0] xx  xx  xx  19
[1] xx  xx  xx  35
[2] xx  xx  xx  72
[3] xx  xx  xx  18

Each line/item represents a person. The first three values are things like income, debt, savings. The fourth item is person age. If you apply decision tree regression with age as the dependent variable, the first split would send items [0], [1], [3] to the root left child, and send item [2] to the root right child, because the ages for items [0], [1], [3] are relatively similar, but the age of 72 for item [2] is much different.

The resulting tree might look like:

                  node0 (root)                  Level 0
             rows [0] [1] [2] [3]

     node1                        node2         Level 1
   rows [0] [2] [3]              row [1]

 node3         node4         node5       node6  Level 2
rows [0] [3]  row [2]        empty       empty

So if age is the dependent variable, a decision tree will place age 72, the anomalous age, in a node that near to the root. Now if you repeat the process, using each column in turn as the dependent variable, and then track which data item appears in the lowest level on average, then you have identified the relatively most anomalous data item.

The scikit-learn library has an IsolationForest module, but I didn’t like several of the details, especially the fact that the scikit module splits randomly instead of using variance reduction. So, I put together my own version that splits in the usual, non-random way. I used C#.

The output of my demo is:

Begin Anomaly Forest demo

Loading synthetic data (200)

First three data items:
 -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

Creating Anomaly Forest object
Done

Analyzing dataset
Using col 0 as dependent variable
Using col 1 as dependent variable
Using col 2 as dependent variable
Using col 3 as dependent variable
Using col 4 as dependent variable
Done

First three anomaly scores:
[0]  8.4000
[1]  10.0000
[2]  9.8000

Most anomalous data item = [122]
Anomaly score = 6.0000

End demo

I called the technique Anomaly Forest to distinguish it from the scikit Isolation Forest.

Ultimately, my Anomaly Forest demo is OK, but because it treats each data column separately, the technique does not take into account interactions between the column variables — so it’s sort of a naive Bayes anomaly detection in a sense. But good fun.

Anomaly detection is like finding a hidden pattern in a dataset. The cover of every issue of Playboy Magazine (except the first one in December 1953), has a bunny logo somewhere. On most covers, the logo is in plain sight, but on some covers, the logo is cleverly hidden.

Left: The logo is disguised as holly leaves on the model’s hat. Right: The logo is disguised as light reflections on the model’s vinyl raincoat, on the lapel.

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 AnomalyForest
{
  internal class AnomalyForestProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin Anomaly Forest demo ");

      // 1. load data
      Console.WriteLine("\nLoading synthetic data (200) ");
      string dataFile =
        "..\\..\\..\\Data\\synthetic_200.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4 };
      double[][] data = MatLoad(dataFile, colsX, ',', "#");
      
      Console.WriteLine("\nFirst three data items: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(data[i], 4, 8);

      Console.WriteLine("\nCreating Anomaly Forest object ");
      AnomalyDetector ad = new AnomalyDetector(10);
      Console.WriteLine("Done ");
      Console.WriteLine("\nAnalyzing dataset ");
      ad.Analyze(data);
      Console.WriteLine("Done ");

      Console.WriteLine("\nFirst three anomaly scores: ");
      for (int i = 0; i "lt" 3; ++i)
        Console.WriteLine("[" + i + "]  " +
          ad.scores[i].ToString("F4"));

      int[] indices = new int[data.Length];
      for (int i = 0; i "lt" indices.Length; ++i)
        indices[i] = i;

      int minIndx = 0; ;
      double minScore = ad.scores[0]; ;
      for (int i = 0; i "lt" ad.scores.Length; ++i)
      {
        if (ad.scores[i] "lt" minScore)
        {
          minScore = ad.scores[i];
          minIndx = i;
        }
      }
      Console.WriteLine("\nMost anomalous data item = [" +
        minIndx + "]");
      Console.WriteLine("Anomaly score = " +
        minScore.ToString("F4"));

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

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

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

    static 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("");
    }

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

  } // class Program

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

  public class AnomalyDetector
  {
    public int maxDepth;  // each tree
    public double[] scores; // one score per data item

    public AnomalyDetector(int maxDepth)
    {
      this.maxDepth = maxDepth;
      this.scores = new double[0];  // null-ish
    }

    public void Analyze(double[][] data)
    {
      int nRows = data.Length;
      int nCols = data[0].Length;
      this.scores = new double[nRows];

      for (int c = 0; c "lt" nCols; ++c)
      {
        Console.WriteLine("Using col " + c +
          " as dependent variable ");
        // make a trainX and trainY
        double[] trainY = new double[nRows];
        double[][] trainX = MatMake(nRows, nCols - 1);
        
        for (int i = 0; i "lt" nRows; ++i)
        {
          int p = 0;  // points into cols of trainX
          for (int j = 0; j "lt" nCols; ++j)
          {
            if (j == c) // at the special column
            {
              trainY[i] = data[i][j];
            }
            else
            {
              trainX[i][p++] = data[i][j];
            }
          } // j
        } // i

        // make a tree minSamples = 2, minLeaf = 1
        DecisionTreeRegressor t = 
          new DecisionTreeRegressor(this.maxDepth, 2, 1, -1);
        t.Train(trainX, trainY);

        // scan tree nodes for assigned rows in leaf nodes
        for (int id = 0; id "lt" t.tree.Count; ++id)
        {
          int level = 
            (int)Math.Truncate(Math.Log2((double)id + 1));
          if (t.tree[id] != null &&
            t.tree[id].rows != null &&
            t.tree[id].isLeaf == true &&
            t.tree[id].rows.Count "gte" 1)
          {
            for (int r = 0; r "lt" t.tree[id].rows.Count; ++r)
            {
              int currRow = t.tree[id].rows[r];
              this.scores[currRow] += level;
            }
          }
        } // each node

      } // c, each column

      // normalize scores relative to number cols
      for (int i = 0; i "lt" this.scores.Length; ++i)
        scores[i] /= nCols;

    } // Analyze()

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

  } // class AnomalyDetector

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

  public class DecisionTreeRegressor
  {
    public int maxDepth;
    public int minSamples;  // aka min_samples_split
    public int minLeaf;  // min number of values in a leaf
    public int numSplitCols; // mostly for random forest
    public List"lt"Node"gt" tree = new List"lt"Node"gt"();
    public Random rnd;  // order in which cols are searched

    public double[][] trainX;  // store data by ref
    public double[] trainY;

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

    public class Node
    {
      public int id;
      public int colIdx;  // aka featureIdx
      public double thresh;
      public int left;  // index into List
      public int right;
      public double value;
      public bool isLeaf;
      public List"lt"int"gt" rows;  // assoc rows in train data

      public Node()
      {
        this.id = -1;
        this.colIdx = -1;
        this.thresh = 0.0;  // aka split value
        this.left = -1;
        this.right = -1;
        this.value = 0.0;  // aka pred y
        this.isLeaf = false;
        this.rows = null;
      }
    } // class Node

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

    public DecisionTreeRegressor(int maxDepth = 2,
      int minSamples = 2, int minLeaf = 1,
      int numSplitCols = -1, int seed = 0)
    {
      // if maxDepth = 0, tree has just a root node
      // if maxDepth = 1, at most 3 nodes (root, l, r)
      // if maxDepth = n, at most 2^(n+1) - 1 nodes
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.numSplitCols = numSplitCols;  // for ran. forest

      // create full tree List with null nodes
      int numNodes = (int)Math.Pow(2, (maxDepth + 1)) - 1;
      for (int i = 0; i "lt" numNodes; ++i)
      {
        this.tree.Add(null);  // empty nodes
      }
      this.rnd = new Random(seed);
    }

    // ------------------------------------------------------
    // public: Train()
    // helpers: MakeTree(), BestSplit(), TreeTargetMean(),
    //   TreeTargetVariance().
    // ------------------------------------------------------

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

    // ------------------------------------------------------
    
    private void MakeTree()
    {
      // no recursion, no pointers, List storage, no stack
      if (this.numSplitCols == -1) // use all cols
        this.numSplitCols = this.trainX[0].Length;

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

      // wait to supply colIdx and thresh in loop
      Node root = new Node();
      root.id = 0;
      root.left = 1;
      root.right = 2;
      root.value = grandMean;
      root.isLeaf = false; // already set
      root.rows = allRows;
      this.tree[0] = root;

      for (int i = 0; i "lt" this.tree.Count; ++i)
      {
        Node currNode = this.tree[i];
        // curr node has values for everything
        // except colIdx and thresh

        // curr node too deep to have children OR
        // curr node not enough rows to split then
        // leave both children as null
        if (currNode == null ||
          currNode.rows.Count == 0) { continue; }

        // if parent cannot be split, make parent a leaf
        if (currNode.id "gte" (int)Math.Pow(2,
          (this.maxDepth)) - 1 ||
          currNode.rows.Count "lt" this.minSamples)
        {
          currNode.isLeaf = true;
          continue;
        }

        // parent has enough rows to try to split
        double[] splitInfo = this.BestSplit(currNode.rows);
        int colIdx = (int)splitInfo[0];
        double splitVal = splitInfo[1]; //split value

        if (colIdx == -1)  // unable split, is a leaf
        {
          currNode.isLeaf = true;
          continue;
        }

        // complete the fields for curr node
        currNode.colIdx = colIdx;
        currNode.thresh = splitVal;

        // construct the children, 
        // except for colIdx and thresh
        // which will be supplied in main loop
        Node leftNode = new Node();
        Node rightNode = new Node();

        // construct children rows using split info
        // all info except colIdx and thresh
        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);
        }

        leftNode.id = currNode.id * 2 + 1;
        if (leftNode.id "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.id = -1;
        leftNode.left = leftNode.id * 2 + 1;
        if (leftNode.left "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.left = -1;
        leftNode.right = leftNode.id * 2 + 2;
        if (leftNode.right "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.right = -1;

        leftNode.rows = leftIdxs;
        leftNode.value =
          this.TreeTargetMean(leftNode.rows);
        this.tree[leftNode.id] = leftNode;

        rightNode.id = currNode.id * 2 + 2;
        if (rightNode.id "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) rightNode.id = -1;
        rightNode.left = rightNode.id * 2 + 1;
        if (rightNode.left "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) rightNode.left = -1;
        rightNode.right = rightNode.id * 2 + 2;
        if (rightNode.right "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) rightNode.right = -1;
        rightNode.rows = rightIdxs;
        rightNode.value =
          this.TreeTargetMean(rightNode.rows);
        this.tree[rightNode.id] = rightNode;

      } // i
      return;
    }

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

    private double[] BestSplit(List"lt"int"gt" rows)
    {
      // implicit params numSplitCols, minLeaf, numSplitCols
      // result[0] = best col idx (as double)
      // result[1] = best split value
      rows.Sort();

      int bestColIdx = -1;  // indicates bad split
      double bestThresh = 0.0;
      double bestVar = double.MaxValue;  // smaller is better

      int nRows = rows.Count;  // or dataY.Length
      int nCols = this.trainX[0].Length;

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

      // process cols in scrambled order
      int[] colIndices = new int[nCols];
      for (int k = 0; k "lt" nCols; ++k)
        colIndices[k] = k;
      // shuffle, inline Fisher-Yates
      int n = colIndices.Length;
      for (int i = 0; i "lt" n; ++i)
      {
        int ri = rnd.Next(i, n);  // be careful
        int tmp = colIndices[i];
        colIndices[i] = colIndices[ri];
        colIndices[ri] = tmp;
      }

      // numSplitCols is usually all columns (-1)
      for (int j = 0; j "lt" this.numSplitCols; ++j)
      {
        int colIdx = colIndices[j];
        HashSet"lt"double"gt" examineds =
          new HashSet"lt"double"gt"();

        for (int i = 0; i "lt" nRows; ++i) // each row
        {
          // if curr thresh been seen, skip it
          double thresh = this.trainX[rows[i]][colIdx];
          if (examineds.Contains(thresh)) continue;
          examineds.Add(thresh);

          // get row idxs where x is lte, gt thresh
          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" nRows; ++k)
          {
            if (this.trainX[rows[k]][colIdx] "lte" thresh)
              leftIdxs.Add(rows[k]);
            else
              rightIdxs.Add(rows[k]);
          }

          // Check if proposed split has too few values
          if (leftIdxs.Count "lt" this.minLeaf ||
            rightIdxs.Count "lt" this.minLeaf)
            continue;  // to next row

          double leftVar =
            this.TreeTargetVariance(leftIdxs);
          double rightVar =
            this.TreeTargetVariance(rightIdxs);
          double weightedVar = (leftIdxs.Count * leftVar +
            rightIdxs.Count * rightVar) / nRows;

          if (weightedVar "lt" bestVar)
          {
            // if this never happens, bestColIdx remains -1
            // which means a bad split. used in MakeTree()
            bestColIdx = colIdx;
            bestThresh = thresh;
            bestVar = weightedVar;
          }

        } // each row
      } // j each col

      double[] result = new double[2];  // out params ugly
      result[0] = 1.0 * bestColIdx;
      result[1] = bestThresh;
      return result;

    } // BestSplit()

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

    private double TreeTargetMean(List"lt"int"gt" rows)
    {
      // mean of rows items in trainY: for node prediction
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
      {
        int r = rows[i];
        sum += this.trainY[r];
      }
      return sum / rows.Count;
    }

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

    private double TreeTargetVariance(List"lt"int"gt" rows)
    {
      double mean = this.TreeTargetMean(rows);
      double sum = 0.0;
      for (int i = 0; i "lt" rows.Count; ++i)
      {
        int r = rows[i];
        sum += (this.trainY[r] - mean) *
          (this.trainY[r] - mean);
      }
      return sum / rows.Count;
    }

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

  } // class DecisionTreeRegressor

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

} // ns

Demo data

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

							

If you have a machine learning regression model that predicts a single numeric value, the three most common ways to evaluate the model are mean squared error (MSE), accuracy, and coefficient of determination (R2).

Note that root mean squared error (RMSE) is just the square root of MSE. This is useful when the target variable to predict has units, such as dollars. MSE has units “dollars-squared” but RMSE has units “dollars”.

MSE is reasonably interpretable, for example, if MSE = 0, the model predicts perfectly. But MSE has no upper limit, and MSE depends on how the data items are scaled.

Accuracy is very interpretable, for example, if accuracy = 75% and there are 200 data items, the model predicts 150 out of 200 correctly. But accuracy requires an arbitrary percentage closeness (perhaps 10%) that determines if a prediction is correct or not.

R2 is sort of a cross between MSE and accuracy. It’s often stated that “R2 is a value between 0 and 1 where higher values indicate a more accurate model.” So, an R2 score of 1 means a model predicts perfectly (which is never possible in practice). R2 is calculated as 1 – (SSres / SStot). The SSres is the sum of the squared differences between target y and predicted y values. The SStot is the sum of squared differences between average of target y and predicted y values. (Note R2 is not the same as classical statistics r2 for correlation).

The R2 score doesn’t isn’t affected by data scaling (good) but R2 can be negative (not good), and there’s no theoretical limit to how negative R2 can be.

As far as I know, all scikit-learn library regression modules (LinearRegression, Ridge, KernelRidge, GradientBoostingRegressor, etc.) have a built in score() function that gives the R2 value for the trained model. Implementing R2 from scratch is easy.

Suppose you have a set of training data. The simplest possible regression model is to just predict the average of the target y values in the training data, for any input x. But if you have a really bad regression model that predicts even worse than just returning the average of the target y values, the R2 score will be negative.

Here’s an example:

x0    x1    y        pred y   ss res   ss tot
0.2   0.3   0.34      0.19    0.0225   0.0841
0.6   0.5   0.86      0.61    0.0625   0.0529
0.3   0.7   0.79      0.44    0.1225   0.0256
0.6   0.2   0.56      0.46    0.0100   0.0049
0.1   0.6   0.61      0.31    0.0900   0.0004
0.3   0.4   0.49      0.29    0.0400   0.0196
0.5   0.5   0.75      0.50    0.0625   0.0144
0.2   0.6   0.64      0.34    0.0900   0.0001
                  
            5.04              0.5000   0.2020
            0.63            
                              R2 = -1.4752   

The target y values are calculated by an unseen y = (x0 * x0) + x1. The predicted y values are computed by the terrible model y’ = (x0 * x0) – x1. R2 = 1 – (0.5000 / 0.2020) = 1 – 2.4752 = -1.4752.

The moral of this blog post is that there’s no single best way to evaluate a regression model.

Evaluating a regression model is mostly objective. Evaluating the name of a hotel is subjective. But here are a couple of names of foreign hotels that don’t resonate well in English.

Demo program.

# ridge_scikit_scratch_r2.py
# scratch R2 with scikit Ridge
# just for fun

import numpy as np
from sklearn.linear_model import Ridge

# Ridge(alpha=1.0, *, fit_intercept=True, copy_X=True,
# max_iter=None, tol=0.0001, solver='auto', positive=False,
# random_state=None)

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

def r2(model, data_X, data_y):
  n = len(data_X)
  ss_res = 0.0
  ss_tot = 0.0
  mean_y = np.mean(data_y)
  for i in range(n):
    x = data_X[i]
    y = data_y[i]
    pred_y = model.predict([x])[0]
    ss_res += (y - pred_y) * (y -pred_y) # inefficient . .
    ss_tot += (y - mean_y) * (y -mean_y)

  r2 = 1.0 - (ss_res / ss_tot) # assume non-zero . .
  return r2

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

print("\nBegin from-scratch R2 demo ")

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

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

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

print("\nCreating scikit Ridge model ")
alpha = 1.0
print("Using default L2 alpha = %0.4f " % alpha)
model = Ridge(alpha=alpha, solver='sag',
  fit_intercept=True, random_state=0)
print("Done ")

print("\nTraining scikit Ridge model ")
print("Using SAG with default training params ")
model.fit(train_X, train_y)
print("Done. Used " + str(model.n_iter_) + " iterations" )

print("\nModel weights: ")
print(model.coef_)
print("Model bias = %0.4f " % model.intercept_)

r2_scikit = model.score(train_X, train_y)
r2_scratch = r2(model, train_X, train_y)

print("\nR2 using built-in score() = %0.4f " % r2_scikit)
print("R2 using from-scratch = %0.4f " % r2_scratch)

print("\nEnd demo ")

Demo 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

Bottom line: I refactored a C# implementation of trimmed kernel ridge regression to approximate support vector regression. The old version had ugly repeated code in two different SGD training phases, so I created a helper method to remove the code repeated-ness.

As is often the case with interesting topics, explaining what the problem is, takes longer than explaining the solution. So bear with me. The goal of a machine learning regression problem is to predict a single numeric value. For example, a bank might want to predict a maximum loan amount based on applicant age, sex, annual income, debt, and so on.

Common regression techniques include linear regression, quadratic regression, nearest neighbors regression, kernel ridge regression, Gaussian process regression, kernel support vector regression, neural network regression, random forest regression, and gradient boost regression. Each technique has dozens of variations, and each technique has pros and cons.

Two closely related machine learning regression techniques are kernel ridge regression (KRR) and support vector regression (SVR). Both techniques use a kernel function (usually the radial basis function, RBF) to compare two data items for similarity. Both techniques must store training data in order to make predictions. KRR must store all training items, while SVR eliminates some of the items during training, leaving just the “support vectors” that need to be stored. So, SVR uses less memory than KRR.

On the other hand, the loss function used by SVR is not differentiable, so SVR cannot be trained using stochastic gradient descent, and therefore, SVR does not scale well to very large datasets. So, KRR is easier to train than SVR.

Note: In addition to the kernel support vector regression technique discussed in this blog post, there is a linear support vector regression technique which is absolutely useless in practice.

One morning before work, I got the idea of combining KRR and SVR. Briefly, I train a KRR model as normal using all data. Then I identify training data items that are predicted “too well” and remove them, leaving just pseudo support vectors. The somewhat counterintuitive SVR idea is that items that are predicted too well don’t help model accuracy very much. Then I retrain a new KRR model using only the reduced training data. This gives a trimmed/sparse KRR model that approximates an SVR model.

In my mind, this technique gives the advantages of KRR (ability to handle very large datasets via SGD training) and the advantages of SVR (fewer stored model items and weights than KRR). I first implemented the idea using the scikit KernelRidge module and the idea worked well. Then I implemented the idea using C#, just to see if I could — using C# I had to implement a lot of code from scratch that is built-in to Python and scikit. The C# code worked but it was clunky so I refactored the clunky C# code to make it more rational. But the refactored C# code was awkward because the SGD part was repeated in the preliminary training phase and the second/final training phase. I refactored the code to avoid the duplicate code.

Let me point out that the idea of trimming a KRR model to approximate an SVR model is fairly obvious, and so I’m reasonably sure other people have used the idea, but I could find no examples or discussion on the Internet.

The output of the new, re-re-factored C# demo is:

Begin kernel ridge regression trim approximation
 to SVR using C#

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 KRR model
Setting RBF gamma = 0.30
Done

Setting lrnRate = 0.0500
Setting maxEpochs = 5000
Setting alpha decay =  0.00001
Setting auto-exit tol = 0.0010
Setting trim epsilon = 0.0035

Training model
epoch =      0 MSE = 0.0181 acc = 0.1700
Auto-exit at epoch 723
epoch =      0 MSE = 0.0377 acc = 0.1009
Auto-exit at epoch 707
Done

Final model weights:
 -0.7406  -0.4260  -0.3825 . . . 
. . .      0.3299   0.8800

Evaluating trimmed KRR model
Number model item/weights = 109

Train acc = 0.9817
Test acc = 0.9500

Train MSE = 0.0001
Test MSE = 0.0003

End KRR trim/prune demo

The demo data is synthetic. There are 200 training items and 40 test items. The demo data was generated by a neural network with random weights and biases. The data is normalized, which is strongly recommended for KRR and SVR because items are compared mathematically and you don’t want a column with huge magnitudes to dominate. Both KRR and SVR can handle non-numeric data, but that’s another story.

The key calling statements are:

double gamma = 0.30;
KRR krr = new KRR(gamma);

double lrnRate = 0.05;
int maxEpochs = 5000;
double alpha = 1.0e-5;  // wt decay regularization
double exitTol = 0.001;  // auto-exit
double epsilon = 0.0035; // trimming
// larger epsilon: fewer wts retained
// smaller epsilon: more wts retained
krr.Train(trainX, trainY, lrnRate, maxEpochs, alpha,
  exitTol, epsilon);

The large number of parameters is a significant weakness of KRR (and SVR too).

After training on 200 data items, those items that were predicted with a small error of less than 0.0035 were removed. There were 91 such well-predicted items, leaving a reduced training dataset of 109 items. These would be called the support vectors in SVR. After training a new KRR model on the reduced dataset, the trimmed/sparse KRR model scored 98.17% accuracy on the reduced training data (107 out of 109 correct) and 95.00% accuracy on the test data (38 out of 40 correct). Very nice.

To validate my C# trimmed KRR approximation to SVR, I ran the data through the scikit SVR module. A grid-search optimized scikit SVR model created 90 support vectors (better than the trimmed KRR 109 items), with 96.50% accuracy on the training data (not quite as good as the C# trimmed KRR model), and 95.00% accuracy on the test data (identical to the trimmed KRR model). In short, the C# trimmed KRR model and the SVR model are very similar. The large number of hyperparameters involved with KRR and SVR makes an exact comparison impossible in practice.

An interesting experiment.

It’s common for me to refactor the software systems I create many times. It’s rarely possible to a get a non-trivial system completely correct on the first effort, and so a system is usually a collection of software sequels, so to speak, where each sequel is a bit better than its predecessor.

Right: In the sequel, “The Chronicles of Riddick” (2004), an alien race, called the Necromongers, are religious fanatics. The story is wildly creative, with very good special effects, great set design, and OK acting. The movie is not well-liked by critics, typically about 1 star out of 4. But I like it a lot and give it my personal B+ grade.

The (re-re-factored) C# trimmed KRR 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 KernelRidgeRegressionTrim
{
  internal class KernelRidgeRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin kernel ridge " +
        "regression trim approximation to SVR using C# ");

      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 KRR model ");
      double gamma = 0.3;    // RBF param
      Console.WriteLine("Setting RBF gamma = " +
        gamma.ToString("F2"));
      KRR krr = new KRR(gamma);
      Console.WriteLine("Done ");

      double lrnRate = 0.05;
      int maxEpochs = 5000;
      double alpha = 1.0e-5;  // wt decay regularization
      double exitTol = 0.001;  // auto-exit
      double epsilon = 0.0035; // trimming
      // larger epsilon: fewer wts retained
      // smaller epsilon: more wts retained
      Console.WriteLine("\nSetting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpochs = " + maxEpochs);
      Console.WriteLine("Setting alpha decay =  " +
        alpha.ToString("F5"));
      Console.WriteLine("Setting auto-exit tol = " +
        exitTol.ToString("F4"));
      Console.WriteLine("Setting trim epsilon = " +
        epsilon.ToString("F4"));
      Console.WriteLine("\nTraining model ");
      krr.Train(trainX, trainY, lrnRate, maxEpochs, alpha,
        exitTol, epsilon);
      Console.WriteLine("Done ");
      Console.WriteLine("\nFinal model weights: ");
      Utils.VecShow(krr.wts, 4, 9);

      Console.WriteLine("\nEvaluating trimmed KRR model ");
      Console.WriteLine("Number model item/weights = " +
        krr.wts.Length);
      double trainAcc = 
        krr.Accuracy(krr.trainX, krr.trainY, 0.10);
      double testAcc = krr.Accuracy(testX, testY, 0.10);

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

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

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

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

  } // class Program

  /*
  Creating scikit SVR model
  Setting gamma = 0.1000
  Setting C = 10.0000
  Setting epsilon = 0.0100

  Number support vectors = 90
  Accuracy (within 0.10) train = 0.9650
  Accuracy (within 0.10) test = 0.9500
  MSE train = 0.0001
  MSE test = 0.0002 
  */

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

  public class KRR
  {
    public double gamma;  // for RBF kernel
    public double[][] trainX;  // need for prediction
    public double[] trainY; // not used this version
    public double[] wts;  // one per trainX item
    public Random rnd; // for SGD

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

    public KRR(double gamma, int seed = 1)
    {
      this.gamma = gamma;
      this.trainX = new double[0][];
      this.trainY = new double[0];
      this.wts = new double[0];
      this.rnd = new Random(seed);  // shuffle train order
    } // ctor

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

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

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

    private void TrainSGD(double[][] trainX, double[] trainY,
     double lrnRate, int maxEpochs, double alpha,
     double exitTol, double epsilon)
    {
      // helper for Train()
      // assign to this.trainX, this.trainY, this.wts

      // trainx, trainY by ref:
      this.trainX = trainX;
      this.trainY = trainY;

      int freq = maxEpochs / 5;  // when to show progress

      // 1. init weights
      this.wts = new double[trainX.Length];
      double lo = -0.10; double hi = 0.10;
      for (int i = 0; i "lt" this.wts.Length; ++i)
        this.wts[i] = (hi - lo) *
          this.rnd.NextDouble() + lo;

      // 2. set up indices for shuffling
      int[] indices = new int[trainX.Length];
      for (int i = 0; i "lt" indices.Length; ++i)
        indices[i] = i;

      // 3. set up prev weights for auto-exit
      double[] prevWeights = new double[trainX.Length];
      for (int j = 0; j "lt" trainX.Length; ++j)
        prevWeights[j] = this.wts[j];

      // 4. SGD prelim train
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        Shuffle(indices);
        for (int i = 0; i "lt" trainX.Length; ++i)
        {
          int idx = indices[i];
          double[] x = trainX[idx];
          double predY = this.Predict(x);
          double actualY = trainY[idx];

          // update wt assoc with x
          this.wts[idx] -= lrnRate * (predY - actualY);
        } // each item

        if (epoch % freq == 0)
        {
          double rmse = this.MSE(trainX, trainY);
          double acc = this.Accuracy(trainX, trainY, 0.10);
          string s1 = "epoch = " +
            epoch.ToString().PadLeft(6);
          string s2 = " MSE = " +
            rmse.ToString("F4");
          string s3 = " acc = " + acc.ToString("F4");
          Console.WriteLine(s1 + s2 + s3);
        }
                
        // if max_change_in_wts / max_weights "lt" tol
        int numWts = this.wts.Length;
        double[] weightDeltas = new double[numWts];
        for (int j = 0; j "lt" numWts; ++j)
          weightDeltas[j] = 
            Math.Abs(prevWeights[j] - this.wts[j]);
        double maxChange = weightDeltas[0];
        for (int j = 0; j "lt" numWts; ++j)
          if (weightDeltas[j] "gt" maxChange)
            maxChange = weightDeltas[j];

        double maxWeightMag = Math.Abs(this.wts[0]);
        for (int j = 0; j "lt" numWts; ++j)
          if (Math.Abs(this.wts[j]) "gt" maxWeightMag)
            maxWeightMag = Math.Abs(this.wts[j]);

        if (maxWeightMag != 0.0 &&
          (maxChange / maxWeightMag) "lt" exitTol)
        {
          Console.WriteLine("Auto-exit at epoch " + epoch);
          break;
        }

        // an auto-exit didn't happen
        for (int j = 0; j "lt" numWts; ++j)
          prevWeights[j] = this.wts[j];
      
      } // each epoch

      // 5. apply one final wt decay regularization
      // as opposed to small decay after every epoch
      for (int j = 0; j "lt" this.wts.Length; ++j)
        this.wts[j] *= (1.0 - alpha);
    }

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

    public void Train(double[][] trainX, double[] trainY,
     double lrnRate, int maxEpochs, double alpha,
     double exitTol, double epsilon)
    {
      // 1. preliminary train
      this.TrainSGD(trainX, trainY, lrnRate,
        maxEpochs, alpha, exitTol, epsilon);

      // 2. construct trimmed data
      int[] isSupportVec = new int[trainX.Length]; 
      for (int i = 0; i "lt" trainX.Length; ++i)
      {
        if (Math.Abs(this.Predict(trainX[i]) - 
          trainY[i]) "gt" epsilon)
          isSupportVec[i] = 1; // just OK prediction
      }
      List"lt"double[]"gt" tmpTrainX = 
        new List"lt"double[]"gt"();
      List"lt"double"gt" tmpTrainY = 
        new List"lt"double"gt"();
      for (int i = 0; i "lt" trainX.Length; ++i)
      {
        if (isSupportVec[i] == 1)
        {
          tmpTrainX.Add(trainX[i]);
          tmpTrainY.Add(trainY[i]);
        }
      }
      double[][] newTrainX = tmpTrainX.ToArray();
      double[] newTrainY = tmpTrainY.ToArray();
      
      // 3. re-train using new small dataset
      this.TrainSGD(newTrainX, newTrainY, lrnRate,
        maxEpochs, alpha, exitTol, epsilon);
    }

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

    private double Rbf(double[] v1, double[] v2)
    {
      // the gamma version aot len_scale version
      int dim = v1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" dim; ++i)
      {
        sum += (v1[i] - v2[i]) * (v1[i] - v2[i]);
      }
      return Math.Exp(-1 * this.gamma * sum);
    }

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

    public double Predict(double[] x)
    {
      int N = this.trainX.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" N; ++i)
      {
        double[] xx = this.trainX[i];
        double k = this.Rbf(x, xx);
        sum += this.wts[i] * k;
      }
      return sum;
    }

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

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

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

    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

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

One of about a dozen common regression techniques is (kernel) support vector regression (SVR). One day before work, I realized that I had not implemented SVR from scratch, using C#, for a long time. So I figured I’d take a stab at it.

There are three main ways to train a support regression model: quadratic programming (QP) optimization, the sequential minimal optimization (SMO) algorithm, and stochastic sub-gradient descent (SSGD / SGD). I used SGD, which is far by the simplest SVR training technique.

The output of my demo is:

Begin C# kernel support vector regression with SGD

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.3000
Setting epsilon = 0.007500
Setting C = 1.00
Setting lrnRate = 0.0010
Setting maxEpochs = 5000
Setting tol = 0.000100
Done

Training SVR model using SGD
epoch =      0 MSE = 0.0430 acc = 0.1300
epoch =   1000 MSE = 0.0001 acc = 0.9850
epoch =   2000 MSE = 0.0001 acc = 0.9850
epoch =   3000 MSE = 0.0001 acc = 0.9800
epoch =   4000 MSE = 0.0001 acc = 0.9800
Done


Model alpha (weights):
 -0.9256  -0.0443  -0.0041  -0.6581  . . .   0.0049
  0.2983   0.4573  -0.0488   0.0956  . . .   0.1471
. . .
  0.0017   0.0023   0.3193   0.0063  . . .  -0.0078
 -0.0690   0.9895   . . .    0.3843
Model bias = 0.4030

Number supp vectors = 194

Evaluating model

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

Train MSE = 0.0001
Test MSE = 0.0001

Train R2 = 0.9985
Test R2 = 0.9949

Predicting for trainX[0]
Predicted y = 0.4922

End SVR with SGD training demo

The demo data is synthetic. It was generated by a 5-10-1 neural network with random weights and bias values. The idea here is that the synthetic data does have an underlying, but complex, non-linear structure which can be predicted.

During training, the SVR model assigns one weight value, into a vector called alpha, to each training item, plus a special weight called the bias. After training, the SVR model determined that 6 of the 200 alpha weights were very close to zero, and so those 6 alpha values and their 6 associated training items were removed. This left 194 weight values and training items, called the support vectors.

All of the parameter values must be determined by trial and error. The gamma parameter defines the RBF function that is used to measure the similarity between data item vectors. Larger values of gamma shrinks the radius of influence of individual training points. This tends to increase model accuracy at the expense of increased risk of model overfitting.

The epsilon value defines how close to correct a prediction must be to be considered a non-support vector. Larger values of epsilon create fewer support vectors.

The C value is used for model regularization, which prevents model alpha weights from becoming very large, which often leads to model overfitting. Larger values of C have a smaller regularization affect.

The lrnRate value controls how much alpha weight values change at each update during training. Larger values of lrnRate increase the speed of training, at the risk of jumping over good weight values.

The maxEpochs value controls how many iterations are performed during training. The effect of larger values of SVR maxEpochs can vary greatly.

The tol (“tolerance”) value controls pruning away training vectors to support vectors, by defining how to close to 0 an alpha weight value must be in order to be pruned away. Larger values of tol allow more alpha weights to be defined as zero, which reduces the number of support vectors.

The biggest weakness of support vector regression is the difficulty of tuning the hyperparameters. Small changes in parameter values can have extremely large changes in the model, and the hyperparameters interact in complex ways.

Support vector regression had a brief surge of popularity in the late 1990s and early 2000s. However, data scientists realized that the closely related kernel ridge regression (KRR) has several significant advantages over SVR, and so the use of SVR declined to the point where it is not used very much today.

SVR is more difficult to implement than KRR, SVR is much more difficult to tune than KRR (KRR can use true SGD, which is easier to tune than SVR sub-gradient descent), and SVR often gives slightly worse prediction accuracy than KRR (due mostly to the difficulty in parameter tuning). That said, there are some problem scenarios where kernel SVR is highly effective.

One of the advantages of implementing a machine learning regression system from scratch is that it gives you a full explanation and understanding of exactly how the system works.

Right: In “Gorath” (1962), a scientist explains that a runaway star named Gorath is on a collision course with Earth. All governments unite to build gigantic rocket thrusters at the South Pole to move Earth out of the way, and then back again. The plan succeeds. My grade for this movie = B.

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

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

// kernel SVR with SGD training
// hard-wired RBF kernel function

namespace SupportVectorRegressionSGD
{
  internal class SupportVectorRegressionSGDProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# kernel support " +
        "vector regression with SGD ");

      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.0075;
      double C = 1.0;
      double lrnRate = 0.001;
      int maxEpochs = 5000;
      double tol = 1.0e-4;

      Console.WriteLine("Setting RBF gamma = " +
        gamma.ToString("F4"));
      Console.WriteLine("Setting epsilon = " +
        epsilon.ToString("F6"));
      Console.WriteLine("Setting C = " +
        C.ToString("F2"));
      Console.WriteLine("Setting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpochs = " +
        maxEpochs);
      Console.WriteLine("Setting tol = " +
        tol.ToString("F6"));

      SVR model = new SVR(gamma, epsilon, C,
        lrnRate, maxEpochs, tol, seed: 0);
      Console.WriteLine("Done ");

      Console.WriteLine("\nTraining SVR model using" +
        " SGD ");
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");
      Console.WriteLine("\nModel alpha (weights): ");
      Utils.VecShow(model.alpha, 4, 9);
      Console.WriteLine("\nModel bias = " + 
        model.b.ToString("F4"));

      Console.WriteLine("\nNumber supp vectors = " +
        model.suppX.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 trainR2 = model.R2(trainX, trainY);
      double testR2 = model.R2(testX, testY);

      Console.WriteLine("\nTrain R2 = " +
        trainR2.ToString("F4"));
      Console.WriteLine("Test R2 = " +
        testR2.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 SGD training demo ");
      Console.ReadLine();
    } // Main()

  } // class Program

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

  public class SVR
  {
    public double gamma;  // for RBF kernel
    public double epsilon;
    public double C; // weight regularization
    public double[][] suppX;  // needed for pred
    public double[] suppY;
    public double[] alpha;  // one per trainX item
    public double b;   // bias
    public double lrnRate;  // for SGD training
    public int maxEpochs;
    public double tol;
    public Random rnd;

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

    public SVR(double gamma, double epsilon, double C,
      double lrnRate, int maxEpochs, double tol,
      int seed = 0)
    {
      this.gamma = gamma;
      this.epsilon = epsilon;
      this.C = C;
      this.suppX = new double[0][]; // compiler happy
      this.suppY = new double[0];
      this.lrnRate = lrnRate;
      this.maxEpochs = maxEpochs;
      this.tol = tol;
      this.alpha = new double[0];
      this.b = 0.0;
      this.rnd = new Random(seed);  // shuffle train order
    } // ctor

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

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

      // init weights
      this.alpha = new double[n];
      double lo = -0.01; double hi = 0.01; // not needed
      for (int i = 0; i "lt" n; ++i)
        this.alpha[i] =
          (hi - lo) * this.rnd.NextDouble() + lo;
      this.b = 0.0;

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

      // set up indices for random order SGD training
      int[] indices = Utils.VecRange(n); // 0, 1, 2, ..
      double lamda = 1.0 / this.C;
      int progressFreq = (int)(this.maxEpochs / 5);

      // main sub-gradient processing loop
      for (int epoch = 0; epoch "lt" this.maxEpochs; ++epoch)
      {
        this.Shuffle(indices);
        for (int i = 0; i "lt" indices.Length; ++i)
        {
          int idx = indices[i];
          double predY = 0.0;
          for (int j = 0; j "lt" this.alpha.Length; ++j)
            predY += this.alpha[j] * K[idx][j]; // fast
          predY += this.b;
          double error = predY - trainY[idx];

          double gradLoss;
          bool insideTube = false;
          if (error "gt" this.epsilon)
            gradLoss = 1.0;
          else if (error "lt" -this.epsilon)
            gradLoss = -1.0;
          else
          {
            gradLoss = 0.0;
            insideTube = true;
          }

          // local kernel regularization gradient
          double gradReg = this.alpha[idx] * K[idx][idx];

          //  decoupled updates to the active index
          this.alpha[idx] -= this.lrnRate *
            (lamda * gradReg + gradLoss);
          this.b -= this.lrnRate * gradLoss;

          // force tiny weights to 0
          if (insideTube == true &&
            Math.Abs(this.alpha[idx]) "lt" this.tol)
            this.alpha[idx] = 0.0;

          // in-loop clip to bound updates mid-flight
          if (this.alpha[idx] "lt" -this.C)
            this.alpha[idx] = -this.C;
          else if (this.alpha[idx] "gt" this.C)
            this.alpha[idx] = this.C;

        } // each item

        // show training progress every few epochs
        if (epoch % progressFreq == 0)
        {
          double mse =
            this.MSE(trainX, trainY);
          double acc =
            this.Accuracy(trainX, trainY, 0.10);
          string s1 = "epoch = " +
            epoch.ToString().PadLeft(6);
          string s2 = " MSE = " +
            mse.ToString("F4");
          string s3 = " acc = " + acc.ToString("F4");
          Console.WriteLine(s1 + s2 + s3);
        }

      } // each epoch

      // final global clip
      for (int i = 0; i "lt" n; ++i)
      {
        if (this.alpha[i] "lt" -this.C)
          this.alpha[i] = -this.C;
        else if (this.alpha[i] "gt" this.C)
          this.alpha[i] = this.C;
      }

      // prune: store only explicit support vectors
      List"lt"int"gt" svLst = new List"lt"int"gt"();
      for (int i = 0; i "lt" this.alpha.Length; ++i)
      {
        if (Math.Abs(this.alpha[i]) "gt" 1.0e-5)
          svLst.Add(i);
      }
      int[] svMask = svLst.ToArray();

      this.suppX = Utils.MatSelectRows(trainX, svMask);
      this.suppY = Utils.VecSelectItems(trainY, svMask);
      this.alpha = Utils.VecSelectItems(this.alpha, svMask);

      return;  // all done
    } // Train

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

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

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

    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(X[i], X[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[] xx = this.suppX[i];
        double k = this.RBF(x, xx);
        sum += this.alpha[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;
    }

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

    public double R2(double[][] dataX, double[] dataY)
    {
      // coefficient of determination
      int n = dataX.Length;
      double ssRes = 0.0; double ssTot = 0.0;
      double meanY = Utils.VecMean(dataY);

      for (int i = 0; i "lt" n; ++i)
      {
        double[] x = dataX[i];
        double actualY = dataY[i];
        double predY = this.Predict(x);
        ssRes += (actualY - predY) * (actualY - predY);
        ssTot += (actualY - meanY) * (actualY - meanY);
      }
      double result = 1.0 - (ssRes / ssTot);
      return result;
    }

  } // 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 n = X.Length;
    //  double[] result = new double[n];
    //  for (int i = 0; i "lt" n; ++i)
    //    result[i] = X[i][col];
    //  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[][] 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 VecMean(double[] vec)
    {
      int n = vec.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += vec[i];
      double result = sum / n;
      return result;
    }

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

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

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

    //public static 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[] 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