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.

Each of the three techniques has several variations, and each variation has several implementation approaches. Briefly, 1.) SGD works for any size of test data but requires you to specify a learning rate and max iterations, which must be determined by trial and error. 2.) left pseudo-inverse works well for small-to-medium size training data but will fail for datasets with unlucky combinations of values (“ill-conditioned”), due to arithmetic overflow or underflow. 3.) MP pseudo-inverse works well for medium-to-large training data, but is more complicated to implement than SGD or left pseudo-inverse.

For MP pseudo-inverse there are six main algorithms (and many secondary algorithms): 1.) SVD (Jacobi algorithm), 2.) SVD (Golub-Kahan algorithm), 3.) SVD (bidiagonalization Householder algorithm), 4.) QR (Givens algorithm), 5.) QR (modified Gram-Schmidt algorithm), 6.) QR (Householder algorithm). If you find this moderately confusing, you are not alone.

I have implemented all of the six MP pseudo-inverse algorithms from scratch. My approach-of-choice is QR-Householder. I recently refactored my QR-Householder implementation to fix a couple of edge-case bugs and remove several inefficiencies. So, I put together a demo of linear regression, trained using my new QR-Householder implementation. Output of the demo:

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

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

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

First three train y:
  0.4840
  0.1568
  0.8054

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

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

Evaluating model

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

MSE train = 0.0026
MSE test = 0.0020

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

Predicted y = 0.5329

End demo

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.

One minor detail. When using MP pseudo-inverse for training, there’s no easy way to introduce regularization. This isn’t a major issue, but if for some reason regularization is absolutely required, you must add a dim-by-dim matrix of zeros, with regularization constant on the diagonal, to the bottom of the training X matrix, and add dim 0.0 values to the end of training y vector. Not worth the trouble in any practical scenario.

I am fascinated to the point of obsession by mathematics, computer science, and machine learning. And for some reason, I love to watch for movie scenes where the closed captioning reads “chittering”. Whenever I find a chittering, I’m thrilled. I’m not joking.

Right: “Harry Potter and the Prisoner of Azkaban” (2004). The third movie in the series is my second favorite. Here Ron’s rat “Scabbers” runs away from Prof. Lupin and Sirius Black. He is revealed to be the evil Peter Pettigrew in disguise.

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

      // 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();
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");

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

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

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

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

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

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

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

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

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

  } // class Program

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

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

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

    public LinearRegressor(int seed = 0)  // ctor
    {
      this.weights = new double[0];
      this.bias = 0;
      this.rnd = new Random(seed); // 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)
    {
      // wts = pinv(designX) * y
      int dim = trainX[0].Length;
      this.weights = new double[dim];

      double[][] X = MatToDesign(trainX);  // design X
      // for regularization, add dim-by-dim with tiny alpha
      // on diagonal to bottom of trainX, and add
      // dim 0.0s to end of trainY
      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 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[M.Length][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols + 1];

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

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

In the back of my mind, I wasn’t entirely happy with my current version of matrix QR decomposition using the Householder algorithm. The current version worked fine, but the underlying code just didn’t feel quite right in some way that I couldn’t articulate.

So I decided to take another stab at matrix decomposition using Householder. For me, the main purpose of QR decomposition is to compute the relaxed Moore-Penrose pseudo-inverse, which in turn is used to train a linear regression model or a quadratic regression model.

After a day or so of work, I got a new version of QR-Householder decomposition up and running to my satisfaction. The output of a demo run:

Begin QR decomposition using Householder algorithm
More precise than QR w/ modified Gram-Schmidt and QR w/
Givens but more complicated.
Latest version (July 2026)

Source (tall) matrix A:
   1.0   2.0   3.0   4.0   5.0
   0.0  -3.0   5.0  -7.0   9.0
   2.0   0.0  -2.0   0.0  -2.0
   4.0  -1.0   5.0   6.0   1.0
   3.0   6.0   8.0   2.0   2.0
   5.0  -2.0   4.0  -4.0   3.0

Computing QR
Done

Q =
  -0.134840   0.258894   0.147094  -0.310903   0.869379
   0.000000  -0.410745   0.759588   0.274531   0.180705
  -0.269680  -0.029872  -0.526675   0.219131   0.300339
  -0.539360  -0.196660   0.116670  -0.738280  -0.260150
  -0.404520   0.776682   0.316260   0.255197  -0.226191
  -0.674200  -0.348511  -0.101841   0.412033   0.049823

R =
  -7.416198  -0.809040  -8.494918  -1.887760  -3.505839
   0.000000   7.303797   2.618812   5.678242  -2.031322
   0.000000   0.000000   7.998637  -2.988836   9.068775
   0.000000   0.000000   0.000000  -8.732743   1.486219
   0.000000   0.000000   0.000000  -0.000000   4.809501

End demo

I validated my QR decomp implementation by sending the same input matrix to the Python language np.linalg.qr() function, to make sure the results were the same.

There are a ton of details. My implementation workks only for matrices that have more rows than columns — such as training data. My implementation returns a ‘reduced’ Q and a ‘reduced’ R (needed for pseudo-inverse) rather than ‘full’ Q and R matrices.

OK. Good fun. Next, I’ll need to implement a relaxed Moore-Penrose pseudo-inverse using my new QR decomp code. And then after that, I’ll need to use the pseudo-inverse to train a linear regression model.

I’m not very good at articulating subjective things. Here are two screen captures from a short AI-generated video titled “Oracle”, from a guy called Anglomangler. I can’t explain why or how, but the nightmarish quality of the video really makes an impact on me.

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

using System;
using System.IO;

namespace MatrixDecompQRHouseholder
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin QR decomposition using" +
        " Householder algorithm ");
      Console.WriteLine("More precise than QR w/ modified" +
        " Gram-Schmidt and QR w/ Givens but more " +
        "complicated. ");
      Console.WriteLine("Latest version (June 2026) ");
      
      double[][] A = new double[6][];
      A[0] = new double[] { 1, 2, 3, 4, 5 };
      A[1] = new double[] { 0, -3, 5, -7, 9 };
      A[2] = new double[] { 2, 0, -2, 0, -2 };
      A[3] = new double[] { 4, -1, 5, 6, 1 };
      A[4] = new double[] { 3, 6, 8, 2, 2 };
      A[5] = new double[] { 5, -2, 4, -4, 3 };

      Console.WriteLine("\nSource (tall) matrix A: ");
      MatShow(A, 1, 6);

      Console.WriteLine("\nComputing QR ");
      double[][] Q;
      double[][] R;
      QRHouseholder.MatDecompQR(A, out Q, out R);
      Console.WriteLine("Done ");

      Console.WriteLine("\nQ = ");
      MatShow(Q, 6, 11);
      Console.WriteLine("\nR = ");
      MatShow(R, 6, 11);

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

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

  } // class Program

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

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

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

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

  } // class QR_Householder

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

} // ns

Python language validation program:

import numpy as np

def main():
  A = np.array([
    [ 1, 2, 3, 4, 5 ],
    [ 0, -3, 5, -7, 9 ],
    [ 2, 0, -2, 0, -2 ],
    [ 4, -1, 5, 6, 1 ],
    [ 3, 6, 8, 2, 2 ],
    [ 5, -2, 4, -4, 3 ]], dtype=np.float64)

  print("\nA = "); print(A)

  # call NumPy function
  print("\nNumPy QR: ")
  Q, R = np.linalg.qr(A, mode='reduced')
  print("\nNumPy Q = "); print(Q)
  print("\nNumPy R = "); print(R)

if name == "__main__":
  main()

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

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

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

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

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

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

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

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

Begin decision tree regression (performant version)

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

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

First three train y:
  0.4840
  0.1568
  0.8054

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

Creating and training tree
Done

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

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

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

MSE train = 0.0048
MSE test = 0.0054

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

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

End demo

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

The diagram below shows how the prediction was arrived at.

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

Good fun.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

      dtr.Explain(x);

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

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

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

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

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

  } // class Program

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        if (currNode == null) break;

        lastValidValue = currNode.value;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  } // class DecisionTreeRegressor

} // ns

Training data:

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

Test data:

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

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

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

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

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

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

Begin prune training data using C# demo

Setting up dummy data

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

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

Prune mask:
  9  5  8

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

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

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

Prune mask:
  7  9  5

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

End demo

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


    //}

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

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

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

  } // class Program

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  } // class Utils

} // ns

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

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

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

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

Scratch linear regression with SGD training

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

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

First three train y:
0.4840
0.1568
0.8054

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

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

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

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

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

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

End demo

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

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

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

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

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

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

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

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

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

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

A fun mental exercise for me on a Sunday morning.

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

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

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

# linear_regression_sgd_regularization_techniques.py

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

import numpy as np

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  print("\nEnd demo ")

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

if __name__ == "__main__":
  main()

Training data:

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

Test data:

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

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

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

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

Here’s the output of my demo:

Begin scratch Python SVR using SMO training

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

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

First three train y:
0.4840
0.1568
0.8054

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

Creating and training SVR model using SMO
Done

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

Model bias = 0.4123

Number support vectors = 197

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

Train MSE = 0.0002
Test MSE = 0.0002

Train R2 = 0.9951
Test R2 = 0.9941

Predicting for train_X[0]
Predicted y = 0.4913

End demo

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

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

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

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

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

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

Good fun.

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

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

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

# svr_smo.py

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

import numpy as np

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  print("\nEnd demo ")

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

if __name__ == "__main__":
  main()

Training data:

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

Test data:

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

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

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

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

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

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

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

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

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

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

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

Begin prune training data demo

Generating dummy training data
Done

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

Pruning down to 3 dissimilar items with RBF gamma = 0.5000

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

End demo

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

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

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

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

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

import numpy as np

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

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

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

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

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

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

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

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

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

  return selected_idxs  

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

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

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

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

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

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

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

  return selected_idxs  

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

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

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

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

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

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

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

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

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

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

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

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

print("\nBegin prune training data demo ")

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

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

gamma = 0.5
n_to_keep = 3

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

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

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

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

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

print("\nEnd demo ")