James D. McCaffrey

I recently (yesterday) made major revisions to my personal C# matrix QR decomposition using the Householder algorithm. In my work environment, I use matrix QR decomposition for matrix pseudo-inverse of a matrix of training data. So, my next step after a new QR decomposition function, presented in yesterday’s blog post, was to implement a new pseudo-inverse function.

Suppose matrix X holds training data (typically a design matrix, where a leading column of 1.0 values has been added to deal with the model bias). Then, to find the pseudo-inverse of X, pinv(X):

1. X = Q * R
2. pinv(X) = inv(Q * R)
3. pinv(X) = inv(Q) * inv(R)
4. pinv(X) = tr(Q) * inv(R)

In step 1, the matrix of training data X is decomposed using QR-Householder decomposition. The QR decomposition applies to any shape matrix, but regular inverse applies only to square matrices. Because training data almost always has more rows than columns, it’s possible to use a relaxed form that assumes this condition.

In step 2, the matrix inverse is applied to both sides of the equation.

In step 3, the matrix property inv(A*B) = inv(B) * inv(A) is applied.

In step 4, the inverse of Q is the transpose of Q, which is very easy to calculate. This is a special property of the Q result of QR decomposition, and is why QR decomposition is used. The R matrix is upper-triangular, and that inverse is a special case which is easy to compute.

Expressed in code:

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

Most of the work is done in the MatDecompQR() function.

Here’s output of a demo run:

Relaxed  MP pseudo-inverse using QR decomposition
 Householder algorithm

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

pseudo-inverse =
   0.0911  -0.0561   0.1650  -0.0250  -0.0144   0.1445
   0.0921  -0.0417   0.0755  -0.1379   0.0939   0.0053
  -0.1618   0.0430  -0.1421   0.1041   0.0790  -0.0414
   0.0664  -0.0250  -0.0145   0.0753  -0.0372  -0.0454
   0.1808   0.0376   0.0624  -0.0541  -0.0470   0.0104

End demo

I partially validated my relaxed MP pseudo-inverse by feeding the same source matrix to the NumPy np.linalg.pinv() function to make sure I got the same result.

OK. Good fun. My next step will be to use this new relaxed MP pseudo-inverse via QR-Householder functionality to train a machine learning regression problem (probably quadratic regression or maybe linear regression).

I don’t collect things. I collect ideas, such as a new algorithm to compute a pseudo-inverse.

Here is card #1.

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 MatrixPseudoInverseQRHouseholder
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nRelaxed  MP pseudo-inverse " +
        "using QR decomposition Householder algorithm ");
            
      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 pinv ");
      double[][] Pinv = QRHouseholder.MatPinv(A);
      Console.WriteLine("Done ");

      Console.WriteLine("\npseudo-inverse = ");
      MatShow(Pinv, 4, 9);
   
      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
  {
    // container class for MP pseudo-inverse via 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

Python validation code:

  import numpy as np

  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)

  Pinv = np.linalg.pinv(A)
  print("\nNumPy pseudo-inverse = "); print(Pinv)

The goal of a machine learning regression problem is to predict a single numeric value. For example, you might want to predict the price of a house in a particular area based on square footage, year built, number of bedrooms, and so on.

There are about a dozen metrics to evaluate a trained regression model, but most are rarely used. In my regression project scenarios, I use prediction mean squared error (MSE), prediction accuracy (to within a specified percentage of the true target value), and coefficient of determination (R2).

A model MSE value is a bit difficult to interpret because 1.) it heavily penalizes outlier predictions, and 2.) it depends on how the target y values are scaled. But MSE is useful because many regression techniques minimize MSE.

Accuracy is the most interpretable metric, but it requires an arbitrary closeness-to-actual-target parameter (like 5% or 10%). Plus accuracy is not very granular, and therefore can be misleading (such as when a model just barely predicts many y values within the specified percentage tolerance).

R2 is sort of a normalized accuracy (larger values are better). R2 doesn’t depend on how the target y values are scaled but R2 isn’t too easy to interpret. R2 is mostly useful to compare two entirely different regression models. R2 is the default “score” metric for scikit-learn regression models and classification models.

A metric called mean absolute scaled error (MASE) is common in the time series regression community, but MASE is rarely used in standard regression problem scenarios. MASE doesn’t heavily penalize outlier predictions like MSE does (because it uses absolute value of error instead of error-squared), and MASE doesn’t depend on how the target y data is scaled.

Note: The MASE implementation shown in this blog post works only for standard regression, not time series regression. The MASE version for TSR predicts a baseline value as the previous target y instead of the average of the target y values.

One Sunday morning, just to entertain myself, I decided to implement a MASE function using the C# language:

  public static double MASE(dynamic model,
    double[][] dataX, double[] dataY)
  {
    // Mean Absolute Scaled Error
    // error relative to baseline predict-mean target y
    // standard tabular data scenarios, not TSR scenarios

    // 1. compute model mean absolute error (MAE)
    int n = dataX.Length;
    double sum = 0.0;
    for (int i = 0; i "lt" n; ++i) // "lt" is less-than
    {
      double predY = model.Predict(dataX[i]);
      sum += Math.Abs(predY - dataY[i]);
    }
    double modelMAE = sum / n;

    // 2. compute mean of target y values
    double sumY = 0.0;
    for (int i = 0; i "lt" n; ++i)
      sumY += dataY[i];
    double meanY = sumY / n;

    // 3. compute baseline MAE (always predict mean y)
    double baseSum = 0.0;
    for (int i = 0; i "lt" n; ++i)
      baseSum += Math.Abs(meanY - dataY[i]);
    double baseMAE = baseSum / n;

    if (baseMAE "lt" 1.0e-12)
      return modelMAE;

    // 4. compute MASE as ratio of model to baseline
    double result = modelMAE / baseMAE;
    return result;
  }

The C# keyword “dynamic” allows the MASE() function to be used with any C# implementation of a regression model (LinearRegression, NearestNeighborsRegressor, QuadraticRegressor, etc.), as long as the model has a Predict() method. For example:

. . .
SomeRegressor model = new SomeRegressor(lrnRate, maxEpochs);
model.Train(trainX, trainY);
maseTrain = MASE(model, trainX, trainY);
Console.WriteLine("MASE for train data = " + maseTrain);
. . .

I put together a demo of MASE with one of my C# linear regression implementations. Output of a demo run:

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

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

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

First three train y:
  0.4840
  0.1568
  0.8054

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

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

Evaluating model

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

MSE train = 0.0026
MSE test = 0.0019

MASE train = 0.2619
MASE test = 0.2402

End demo

Because MASE is an error metric, for MASE, lower values are better. The MASE metric can be as low as 0.0 (perfect predictions, no error — impossible in practice) up to 1.0 (the model predicts exactly like just predicting the average of the target y values for any input), or MASE can be greater than 1.0 which means the model predicts even worse than just predicting the average y for any input.

MASE has some nice properties. And MASE is very common in the time series regression community, but MASE is almost never used in the tabular data regression community. The main reason why R2 is nearly universal for tabular data regression is mostly because R2 is used by scikit-learn, and scikit-learn has a virtual monopoly on machine learning library code and usage patterns.

See https://jamesmccaffreyblog.com/2026/04/20/refactoring-my-demo-of-linear-regression-with-closed-form-training-using-csharp/ for one version of linear regression code, and the data.

I am a big fan of old science fiction movies from the 1950s and 1960s. The limited non-digital technology of the time meant that all movies were relatively expensive and difficult to create, and so there aren’t many sci-fi movies from those two decades. But starting in the 1990s, digital technology advanced to the point where a reasonably professional movie could be made for roughly $2.0 million dollars. This led to an explosion of sci-fi movies. Most of these movies are bad, but every now and then a low-budget sci-fi movie surprises me in a good way. Unlike a regression model where evaluation is objective, it’s difficult to compute a subjective evaluation metric for movies.

In “Ice Spiders” (2007), a group of young skiers on an isolated Utah mountain run into giant genetically altered spiders that are the product of a secret government lab. This movie has a nice combination of interesting plot, pretty good acting, humor-vs-scary, and excellent special effects. The movie doesn’t try to be anything more than it is. My grade = solid B.

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