James D. McCaffrey

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

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

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

Here’s the output of my demo:

Begin scratch Python SVR using SMO training

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

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

First three train y:
0.4840
0.1568
0.8054

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

Creating and training SVR model using SMO
Done

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

Model bias = 0.4123

Number support vectors = 197

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

Train MSE = 0.0002
Test MSE = 0.0002

Train R2 = 0.9951
Test R2 = 0.9941

Predicting for train_X[0]
Predicted y = 0.4913

End demo

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

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

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

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

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

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

Good fun.

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

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

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

# svr_smo.py

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

import numpy as np

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  print("\nEnd demo ")

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

if __name__ == "__main__":
  main()

Training data:

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

Test data:

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

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

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

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

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

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

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

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

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

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

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

Begin prune training data demo

Generating dummy training data
Done

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

Pruning down to 3 dissimilar items with RBF gamma = 0.5000

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

End demo

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

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

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

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

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

import numpy as np

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

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

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

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

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

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

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

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

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

  return selected_idxs  

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

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

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

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

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

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

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

  return selected_idxs  

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

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

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

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

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

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

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

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

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

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

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

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

print("\nBegin prune training data demo ")

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

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

gamma = 0.5
n_to_keep = 3

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

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

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

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

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

print("\nEnd demo ")

An Extra Trees Regressor (“Extremely Randomized Trees”) is an ensemble machine learning algorithm that builds multiple randomized decision trees to predict continuous numerical values. It differs from Random Forest regressions by introducing greater randomization during node splitting and utilizing the entire original dataset to train each tree rather than using samples of the original dataset.

One weekend I figured I’d implement Extra Trees Regression, from scratch, using Python and NumPy. The process took me a few hours, but was manageable because I have quite a bit of experience with tree-based systems.

The output of my demo:

Extra trees regression scratch Python NumPy

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

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

First three y targets:
0.4840
0.1568
0.8054

Setting n_estimators = 200
Setting max_depth = 6
Setting max_features = 4

Accuracy on train (within 0.10) = 0.7700
Accuracy on test (within 0.10) = 0.7250

MSE on train = 0.0006
MSE on test = 0.0012

Predicting for:
[[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]]
Predicted y = 0.4930

SCIKIT:

Accuracy on train (within 0.10) = 0.7750
Accuracy on test (within 0.10) = 0.6750

MSE on train = 0.0006
MSE on test = 0.0012

Done

The data is synthetic. It was generated by a neural network with random weights and biases. There are five predictor values and one value to predict. There are 200 training items and 40 test items.

Extra Trees Regression is an ensemble technique. For my demo, I used a collection of 200 base decision trees. The final prediction is the average of the 200 base trees. The base trees are exactly like standard decision trees, except, during construction, at each node, a base Extra Tree picks one value at random from each active column, and then uses the one value that gives the best split — a separation of rows of data that has the lowest variance (so that the left child nodes and right child nodes have relatively similar target y values).

The ideas are suprisingly subtle and implementation is extremely tricky.

The output of my from-scratch demo is very close to the output from the scikit-learn ExtraTreesRegressor module. The output wasn’t expected to be identical because Extra Trees Regression is fundamentally random.

Extra Trees Regression has no super powers — it’s just math and computer science. When I was a young man, I loved, and learned to read from, comic books of the early 1960s. I especially loved heros who had no special powers — they relied on just their intelligence — Batman, Green Arrow, Adam Strange, and Iron Man.

Right: In issue #48 (December 1963), the suit changed to red and gold. Iron Man got his own title in May 1968.

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

# extra_trees_regression_scratch.py

import numpy as np

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

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

class MyExtraTreesRegressor:
  def __init__(self, n_estimators=100, max_depth=3,
    min_samples_split=2, max_features='auto',
    RandomState=0):

    self.n_estimators = n_estimators
    self.max_depth = max_depth
    self.min_samples_split = min_samples_split
    self.max_features = max_features
    self.RandomState = RandomState
    self.trees = []
    self.rng = np.random.RandomState(RandomState)

  def fit(self, train_X, train_y):
    X = np.asarray(train_X)
    n_samples, n_features = X.shape

    # map scikit-learn standard max_features to your
    # n_split_cols parameter
    if self.max_features == 'auto' or self.max_features is None:
      n_split_cols = n_features
    elif self.max_features == 'sqrt':
      n_split_cols = int(np.sqrt(n_features))
    elif self.max_features == 'log2':
      n_split_cols = int(np.log2(n_features))
    else:
      n_split_cols = int(self.max_features)
    n_split_cols = max(1, min(n_split_cols, n_features))

    self.trees = []
    for _ in range(self.n_estimators):
      tree_seed = self.rng.randint(0, 1000000)
      
      # instantiate our specialized base tree 
      tree = MyExtraTreesBaseTree(
        max_depth=self.max_depth,
        min_samples=self.min_samples_split,
        min_leaf=1,
        n_split_cols=n_split_cols,
        save_rows=False,
        seed=tree_seed
      )
      tree.fit(train_X, train_y)
      self.trees.append(tree)
    return self

  def predict(self, X):
    X = np.asarray(X)
    preds = np.array([tree.predict(X) for tree in self.trees])
    return np.mean(preds, axis=0)


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

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

def accuracy(model, data_X, data_y, pct_close):
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    y_pred = model.predict(x)[0]

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

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

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

  return sum / n

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

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

class MyExtraTreesBaseTree:

  # identical to MyDecisionTreeRegressor except best_split()

  def __init__(self, max_depth=3, min_samples=2, 
    min_leaf=1, n_split_cols=-1, save_rows=False, 
    seed=0):

    self.max_depth = max_depth
    self.min_samples = min_samples   # to split
    self.min_leaf = min_leaf         # after split
    self.n_split_cols = n_split_cols # -1 all cols
    self.save_rows = save_rows
    self.rnd = np.random.RandomState(seed)

    self.tree = []  # list: no recursive, no ptrs
    num_nodes = 2**(max_depth + 1) - 1
    for i in range(num_nodes):
      self.tree.append(None)

    self.train_X = None
    self.train_y = None

  # ...............................................

  class Node:
    def __init__(self):
      self.id = -1
      self.col_idx = -1  # aka split col
      self.thresh = 0.0  # aka split val
      self.left = -1
      self.right = -1
      self.value = 0.0   # predicted y
      self.is_leaf = False
      self.rows = []     # associated rows in train

  # ...............................................

  def fit(self, train_X, train_y):
    # aka train(), build_tree()
    # self.train_X = np.array(train_X)  # safety
    # self.train_y = np.array(train_y)
    self.train_X = train_X
    self.train_y = train_y
    
    max_id = 2**(self.max_depth+1) - 2  # md = 3, m_id = 14
    max_start_id = 2**self.max_depth - 1  # md = 3, ms = 7

    # prep root node
    all_rows = []
    for i in range(len(self.train_X)):  # or use arange
      all_rows.append(i)
    grand_mean = self.tree_target_mean(all_rows)
    
    root = self.Node()
    root.id = 0
    root.value = grand_mean
    root.is_leaf = False
    root.rows = all_rows
    self.tree[0] = root

    for i in range(len(self.tree)):  # each node
      curr_node = self.tree[i]       # convenience
    
      if curr_node is None: continue

      # if node too deep to have children
      #  or not enough rows to split, leave children alone
      if curr_node.id "gte" max_start_id or \
        len(curr_node.rows) "lt" self.min_samples:
        curr_node.is_leaf = True
        continue

      # try to split curr node
      col_idx, split_val = self.best_split(curr_node.rows)
      if col_idx == -1:  # bad split
        curr_node.is_leaf = True
        curr_node.left = -1
        curr_node.right = -1
        continue
  
      # got good split info
      curr_node.col_idx = col_idx
      curr_node.thresh = split_val

      # make rows for the children
      left_idxs = []
      right_idxs = []
      for k in range(len(curr_node.rows)):
        r = curr_node.rows[k]
        if self.train_X[r][col_idx] "lte" split_val:
          left_idxs.append(r)
        else:
          right_idxs.append(r)

      # make left child
      left_id = curr_node.id * 2 + 1
      if left_id "lte" max_id and len(left_idxs) "gte" \
        self.min_leaf:
        curr_node.left = left_id
        left_node = self.Node()
        left_node.id = left_id
        left_node.rows = left_idxs
        left_node.value = \
          self.tree_target_mean(left_node.rows)
        self.tree[left_id] = left_node
      else:
        curr_node.left = -1

      # make right child
      right_id = curr_node.id * 2 + 2
      if right_id "lte" max_id and len(right_idxs) "gte" \
        self.min_leaf:
        curr_node.right = right_id
        right_node = self.Node()
        right_node.id = right_id
        right_node.rows = right_idxs
        right_node.value = \
          self.tree_target_mean(right_node.rows)
        self.tree[right_id] = right_node
      else:
        curr_node.right = -1

      if curr_node.left == -1 and curr_node.right == -1:
        curr_node.is_leaf = True

    # zap away rows to save space if in ensemble
    if self.save_rows == False:
      for k in range(len(self.tree)):
        if self.tree[k] is not None:
          self.tree[k].rows = None

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

  def predict_one(self, x):
    # x is a vector
    curr_idx = 0
    last_valid_value = 0.0
        
    while curr_idx != -1 and curr_idx "lt" len(self.tree):
      curr_node = self.tree[curr_idx]
      
      if curr_node is None: break  # safety check
                
      last_valid_value = curr_node.value
      if curr_node.is_leaf == True: break

      if x[curr_node.col_idx] "lte" curr_node.thresh:
        curr_idx = curr_node.left
      else:
        curr_idx = curr_node.right

    return last_valid_value

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

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

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

  def best_split(self, rows):  # Extra Trees version
    best_col_idx = -1
    best_thresh = 0.0
    best_var = float('inf')
    n_rows = len(rows)
    n_cols = len(self.train_X[0])

    # Fisher-Yates shuffle to get candidate features
    col_indices = np.arange(n_cols)
    for i in range(n_cols - 1):
      ri = self.rnd.randint(i, n_cols)
      tmp = col_indices[i]
      col_indices[i] = col_indices[ri]
      col_indices[ri] = tmp

    n_cols_to_use = min(self.n_split_cols, n_cols) \
      if self.n_split_cols != -1 else n_cols
    active_cols = col_indices[0:n_cols_to_use]

    # evaluate a threshold for each active column
    for col_idx in active_cols:
      feature_vals = self.train_X[rows, col_idx]
      min_v = np.min(feature_vals)
      max_v = np.max(feature_vals)
      
      # if all feature values are identical, 
      # this column cannot be split
      if min_v == max_v:
        continue
      
      # EXTRA TREES MECHANIC: 
      # select a single threshold completely at random
      thresh = self.rnd.uniform(min_v, max_v)
      
      left_count = 0
      right_count = 0
      left_sum_y = 0.0
      left_sum_sq_y = 0.0
      
      # count and sum target values to compute
      # variance reduction
      for r in rows:
        y_curr = self.train_y[r]
        if self.train_X[r, col_idx] "lte" thresh:
          left_count += 1
          left_sum_y += y_curr
          left_sum_sq_y += y_curr * y_curr
        else:
          right_count += 1

      # enforce min_leaf boundaries
      if left_count "lt" self.min_leaf or \
        right_count "lt" self.min_leaf:
        continue

      # isolate right side aggregates dynamically
      # using total sums
      total_sum_y = np.sum(self.train_y[rows])
      total_sum_sq_y = np.sum(self.train_y[rows] ** 2)
      
      right_sum_y = total_sum_y - left_sum_y
      right_sum_sq_y = total_sum_sq_y - left_sum_sq_y

      # variance shortcut formula: E[X^2] - (E[X])^2
      left_var = (left_sum_sq_y / left_count) - \
        ((left_sum_y / left_count) ** 2)
      right_var = (right_sum_sq_y / right_count) - \
        ((right_sum_y / right_count) ** 2)
      
      if left_var "lt" 0.0: left_var = 0.0
      if right_var "lt" 0.0: right_var = 0.0
      
      weighted_var = ((left_count * left_var) + \
        (right_count * right_var)) / n_rows
      
      if weighted_var "lt" best_var:
        best_var = weighted_var
        best_col_idx = col_idx
        best_thresh = thresh

    return best_col_idx, best_thresh

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

  def tree_target_mean(self, rows):
    if rows is None or len(rows) == 0: return 0.0
    sum = 0.0
    for i in range(len(rows)):
      r = rows[i]
      sum += self.train_y[r]
    return sum / len(rows)

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

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

print("\nExtra trees regression scratch Python NumPy ")

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

test_file = ".\\Data\\synthetic_test_40.txt"
test_X = np.loadtxt(test_file, comments="#",
  usecols=[0,1,2,3,4],
  delimiter=",",  dtype=np.float64)
test_y = np.loadtxt(test_file, comments="#", usecols=5,
  delimiter=",",  dtype=np.float64)
print("Done ")

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

n_ests = 200
max_d = 6
max_f = 4

print("\nSetting n_estimators = " + str(n_ests))
print("Setting max_depth = " + str(max_d))
print("Setting max_features = " + str(max_f))

model = MyExtraTreesRegressor(n_estimators=n_ests,
 max_depth=max_d, min_samples_split=2, 
 max_features=max_f, RandomState=0)

model.fit(train_X, train_y)

acc_train = accuracy(model, train_X, train_y, 0.10)
print("\nAccuracy on train (within 0.10) = \
%0.4f " % acc_train)
acc_test = accuracy(model, test_X, test_y, 0.10)
print("Accuracy on test (within 0.10) = \
%0.4f " % acc_test)

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

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

# ==== SCIKIT ===========

from sklearn.ensemble import ExtraTreesRegressor

print("\nSCIKIT: ")

model = ExtraTreesRegressor(n_estimators=n_ests, 
  max_depth=max_d, max_features=max_f, random_state=0)
model.fit(train_X, train_y)

acc_train = accuracy(model, train_X, train_y, 0.10)
print("\nAccuracy on train (within 0.10) = \
%0.4f " % acc_train)
acc_test = accuracy(model, test_X, test_y, 0.10)
print("Accuracy on test (within 0.10) = \
%0.4f " % acc_test)

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

print("\nDone ")

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 machine learning, if training data is multicollinear, the interpretability of the resulting model will likely be poor (the details are tricky and out of the scope of this post). Multicollinear data has two or more rows that are highly correlated mathematically. For example, if a set of training data has a column (predictor variable) house-size (in square feet), and another column house-price, those two columns are likely highly correlated with each other.

Note: Multicollinear training data is usually not a problem for neural network regression, and not a problem for tree-based regression (random forest, gradient boosting).

The most common way to analyze training data for multicollinearity is to compute the VIF (variance inflation factor) for each column of the data.

VIF is a value between 1.0 and positive infinity (in weird scenarios, a VIF value could be less than 1.0). Briefly, if all column VIF values are less than about 7.0, the data is probably OK in terms of multicollinearity.

if VIF is close to 1.0, the column is not correlated with other cols
if VIF between 1.0 and 5.0, column is mildly correlated
if VIF between 5.0 and 10.0, column is highly correlated
if VIF greater than 10.0, column is extremely correlated

To compute the VIF for a specified column of training data, you use the specified column as the dependent y variable, and use the remaining columns as the independent predictor variables, and compute a linear regression model, and then compute the R2 (coefficient of determination) for the model. The VIF value for the column is 1.0 / (1.0 – R2).

Suppose that you have a set of training data X predictor values, and you use some column c as the dependent y variable, and all the other columns as predictors for c. After training a linear regression model, you compute R2 and it is 0.90 — which means column c is predicted very well by the other columns. The VIF value for column c is 1.0 / (1.0 – R2) = 1.0 / 0.10 = 10.0 which is large which is bad because column c is a linear combination of the other columns — the data is somewhat multicollinear.

Now, with the same setup, suppose R2 is 0.20 — which means column c cannot be predicted well by the other columns. The VIF value is 1.0 / (1.0 – 0.20) = 1.0 / 0.8 = 1.25 which is a small value, which is good, because column c is not a linear combination of the other columns, and therefore the data is not multicollinear.

Some time ago, I put together a demo using Python NumPy and the scikit library. This was easy because NumPy and scikit had all the components available.

I decided to entertain myself by refactoring my Python demo to the C# language. This required a bit of effort because I had to implement a no-touch linear regression class, plus several helper functions, such as a program-defined MatWithoutColumn() to mimic the built-in numpy np.delete() function.

I created two datasets. The first data set has five columns of predictors, followed by a column of target y values. The data is “normal” in the sense that there’s no multicollinearity. There are 20 items. 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 second dataset is highly multicollinear, where the third column is 2 times the first column, plus the second column, plus a small random value between 0.000 and 0.001. It looks like:

-0.1660,  0.4406,  0.1096, -0.3953, -0.7065, 0.4840
 0.0776, -0.1616, -0.0045, -0.5911,  0.7562, 0.1568
-0.9452,  0.3409, -1.5482,  0.1174, -0.7192, 0.8054
. . .

The output of my C# VIF demo program is:

Begin variance inflation factor (VIF) demo using C#

Loading synthetic (20) normal data

First three lines:
  -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

Begin VIF analysis
col =  0  |  vif = 1.1979
col =  1  |  vif = 1.4590
col =  2  |  vif = 1.2345
col =  3  |  vif = 1.3025
col =  4  |  vif = 1.2120

Loading synthetic (20) highly multicollinear data
(col[2] = 2.0 * col[0] + col[1] + noise)

First three lines:
  -0.1660   0.4406   0.1096  -0.3953  -0.7065
   0.0776  -0.1616  -0.0045  -0.5911   0.7562
  -0.9452   0.3409  -1.5482   0.1174  -0.7192

Begin VIF analysis
col =  0  |  vif = 25101680.1917
col =  1  |  vif = 5769710.7907
col =  2  |  vif = 30816932.6001
col =  3  |  vif = 1.2937
col =  4  |  vif = 1.1879

End demo

As expected, the first dataset didn’t have any bad VIF values. But the VIF values for the second dataset show that columns [0], [1], [2] are highly correlated.

For my internal linear regressor, I used SGD training with some default parameters and auto-early-exit, so that the regressor just works without tuning. While experimenting, I noticed that SGD training has a lot of trouble dealing with multicollinear data — which is the whole point of checking for multicollinear data in the first place.

An interesting exploration.

In machine learning, you don’t want a relationship between two columns in your training data. But in science fiction movies, you absolutely do want a good relationship between the hero and the main actress.

Two huge lost opportunities because of no actor personality multicollinearity.

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

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

namespace VarianceInflationFactor
{
  internal class VarianceInflationFactorProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin variance inflation " +
        "factor (VIF) demo using C# ");

      Console.WriteLine("\nLoading synthetic (20)" +
        " normal data ");
      string file1 = 
        "..\\..\\..\\Data\\synthetic_train_20.txt";
      double[][] train1 = MatLoad(file1, 
        new int[] { 0, 1, 2, 3, 4 }, ',', "#");
      Console.WriteLine("\nFirst three lines: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(train1[i], 4, 9);

      Console.WriteLine("\nBegin VIF analysis ");
      for (int j = 0; j "lt" train1[0].Length; ++j)
      {
        double z = VarInfFactor(train1, j);
        Console.WriteLine("col = " +
          j.ToString().PadLeft(2) + 
          "  |  vif = " + z.ToString("F4"));
      }

      Console.WriteLine("\nLoading synthetic (20)" +
        " highly multicollinear data ");
      Console.WriteLine("(col[2] = 2.0 * col[0] +" +
        " col[1] + noise) ");
      string file2 = 
        "..\\..\\..\\Data\\synthetic_train_20_collinear.txt";
      double[][] train2 = MatLoad(file2,
        new int[] { 0, 1, 2, 3, 4 }, ',', "#");
      Console.WriteLine("\nFirst three lines: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(train2[i], 4, 9);

      Console.WriteLine("\nBegin VIF analysis ");
      for (int j = 0; j "lt" train2[0].Length; ++j)
      {
        double z = VarInfFactor(train2, j);
        Console.WriteLine("col = " +
          j.ToString().PadLeft(2) +
          "  |  vif = " + z.ToString("F4"));
      }
  
      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

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

    static double VarInfFactor(double[][] data, int col)
    {
      // predict col in data[][] using other columns
      double[][] X = MatWithoutCol(data, col);
      double[] y = MatGetColumn(data, col);
      LinearRegressor model = new LinearRegressor();
      model.Train(X, y);
      //model.TrainLeftPinv(X, y);
      double r2 = model.R2(X, y);  // usually in (0.0, 1.0)
      double vif = 1.0 / (1.0 - r2);  // r2 could be 1
      return vif;
    }

    // ------------------------------------------------------
    // helpers for Main(): MatLoad(), 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 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("");
    }

    // ------------------------------------------------------
    // helpers for VarInflationFactor()
    // ------------------------------------------------------

    static double[][] MatWithoutCol(double[][] data, int col)
    {
      int nRows = data.Length;
      int nCols = data[0].Length;

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

      int k = 0; // into result
      for (int j = 0; j "lt" nCols; ++j)
      {
        if (j == col) continue;
        for (int i = 0; i "lt" nRows; ++i)
          result[i][k] = data[i][j];
        ++k;
      }
      return result;
    }

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

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

  } // class Program

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

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

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

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

    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 int Train(double[][] trainX, double[] trainY,
      double initRate=0.10, int maxEpochs=10000,
      double noChangeTol=0.00001, int consecutiveNoChange=5)
    {
      // exit when dist(old wts, new wts) lt noChangeTol
      // for consecutiveNoChange times
      int n = trainX.Length;
      int dim = trainX[0].Length;
      this.weights = new double[dim];

      // initialize weights and bias
      // technically not necessary
      double low = -0.01; double hi = 0.01;
      for (int i = 0; i "lt" dim; ++i)
        this.weights[i] = (hi - low) *
          this.rnd.NextDouble() + low;
      this.bias = (hi - low) *
          this.rnd.NextDouble() + low;

      int[] indices = new int[n];  // of train data
      for (int i = 0; i "lt" n; ++i)
        indices[i] = i;

      double[] oldWeightsAndB = new double[dim+1];
      double[] newWeightsAndB = new double[dim+1];
      int countNoChange = 0;
      
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        Shuffle(indices, this.rnd);

        double lrnRate = 
          initRate / Math.Pow((double)(epoch+1), 0.25);
        // Console.WriteLine(lrnRate);

        for (int j = 0; j "lt" dim; ++j)
          oldWeightsAndB[j] = this.weights[j];
        oldWeightsAndB[dim] = this.bias;

        for (int i = 0; i "lt" n; ++i) // each train item
        {
          int idx = indices[i];
          double[] x = trainX[idx];
          double predY = this.Predict(x);
          double actualY = trainY[idx];

          // decay weights before update
          // not needed because regularization not useful
          // for (int j = 0; j "lt" dim; ++j) // each weight
          //   this.weights[j] *= (1 - decay);

          // update weights and bias
          for (int j = 0; j "lt" dim; ++j) // each weight
            this.weights[j] -= lrnRate *
              (predY - actualY) * x[j];
          this.bias -= lrnRate * (predY - actualY) * 1;
        }

        //// display progress 5 times: very noisy
        //if (epoch % (int)(maxEpochs / 5) == 0) // progress
        //{
        //  double r2 = this.R2(trainX, trainY);
        //  string s = "";
        //  s += "epoch = " + epoch.ToString().PadLeft(5);
        //  s += "  R2 = " + r2.ToString("F4").PadLeft(8);
        //  Console.WriteLine(s);
        //}

        // check for early-exit after each epoch
        for (int j = 0; j "lt" dim; ++j)
          newWeightsAndB[j] = this.weights[j];
        newWeightsAndB[dim] = this.bias;
        double d = EuclideanDist(oldWeightsAndB,
          newWeightsAndB);

        if (d "lt" noChangeTol)
        {
          ++countNoChange;
          if (countNoChange == consecutiveNoChange)
            return epoch;
        }
        else
         countNoChange = 0; // reset

      } // epoch
      return maxEpochs;

    } // Train

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

    private static double EuclideanDist(double[] v1,
      double[] v2)
    {
      int n = v1.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += (v1[i] - v2[i]) * (v1[i] - v2[i]);
      return Math.Sqrt(sum);
    }

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

    public double R2(double[][] dataX, double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;

      for (int i = 0; i "lt" n; ++i)
        sum += dataY[i];
      double meanY = sum / n;

      double ssRes = 0.0; // sum squared residuals
      double ssTot = 0.0; // sum squared total
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes +=
          (dataY[i] - predY) * (dataY[i] - predY);
        ssTot +=
          (dataY[i] - meanY) * (dataY[i] - meanY);
      }
      return 1.0 - (ssRes / ssTot);
    }

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

    private static void Shuffle(int[] indices, Random rnd)
    {
      int n = indices.Length;
      for (int i = 0; i "lt" n; ++i) // one extra pass
      {
        int ri = rnd.Next(i, n);
        int tmp = indices[i];
        indices[i] = indices[ri];
        indices[ri] = tmp;
      }
    }

  } // class LinearRegressor

} // ns

First, normal, no multicollinearity dataset:

# synthetic_train_20.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

Second, highly multicollinear, dataset:

# synthetic_train_20_collinear.txt
# col [2] = 2*[0] + [1] + rand(0.001)
#
-0.1660,  0.4406,  0.1096, -0.3953, -0.7065, 0.4840
 0.0776, -0.1616, -0.0045, -0.5911,  0.7562, 0.1568
-0.9452,  0.3409, -1.5482,  0.1174, -0.7192, 0.8054
 0.9365, -0.3732,  1.5016,  0.7528,  0.7892, 0.1345
-0.8299, -0.9219, -2.5800,  0.7563, -0.8033, 0.7955
 0.0663,  0.3838,  0.5179,  0.3730,  0.6693, 0.3206
-0.9634,  0.5003, -1.4245,  0.4963, -0.4391, 0.7377
-0.1042,  0.8172,  0.6100, -0.4244, -0.7399, 0.4801
-0.9613,  0.3577, -1.5636, -0.4689, -0.0169, 0.6861
-0.7065,  0.1786, -1.2325, -0.7953, -0.1719, 0.5569
 0.3888, -0.1716,  0.6073,  0.0718,  0.3276, 0.2500
 0.1731,  0.8068,  1.1544, -0.7214,  0.6148, 0.3297
-0.2046, -0.6693, -1.0770, -0.3045,  0.5016, 0.2129
 0.2473,  0.5019,  0.9980, -0.4601,  0.7918, 0.2613
-0.1438,  0.9297,  0.6435,  0.2434, -0.7705, 0.5171
 0.1568, -0.1837,  0.1313,  0.8068,  0.1474, 0.3307
-0.9943,  0.2343, -1.7528,  0.0541,  0.7719, 0.5581
 0.2467, -0.9684, -0.4732,  0.3818,  0.9946, 0.1092
-0.6553, -0.7257, -2.0345,  0.3936, -0.8680, 0.7018
 0.8460,  0.4230,  2.1166, -0.9602, -0.9476, 0.1996

I set out to implement VIF (variance inflation factor) from scratch using C#. VIF is a metric that indicates if there is multicollinearity in a set of training data. If there are n columns of predictors, then you need to compute n linear regression models and their R2 scores.

But the linear regression models just have to work without any hyperparameter tuning. Therefore, I needed to implement such a no-tune linear regression system using C#. I decided to use SGD training (stochastic gradient descent) with an adaptive learning rate, and automatic early-exit logic. Note: For small to medium datasets, it’s probably better to use a closed-form training algorithm, specifically Moore-Penrose pseudo-inverse.

There are many ways to implement both ideas. For the adaptive learning rate, I used the scheme in the scikit SGDRegressor module. An initial, somewhat large, learning rate of 0.01 is set. Then, at each training epoch t, lr(t) = 0.01 / t^0.25. This will slowly reduce the learning rate.

For the auto-exit from training, I used a scheme that I came across but I don’t remember when/where. The stopping condition is when (max_change_in_wts / max_weight) less-than 0.001. This is simple but OK, because the linear regression models for VIF don’t need to be optimal in terms of minimizing mean squared error. Note: Another possibility is to track change in Euclidean distance of a vector holding weights + bias.

For my no-tune SGD linear regression demo, I used one of my standard synthetic 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
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
. . .

The data was generated by a 5-10-1 neural network with random weights and biases. There are 200 training items and 40 test items.

The output of my demo is:

Begin C# linear regression SGD training with
 adaptive LR and auto-exit

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
Done
Used 179 epochs

Weights/coefficients:
-0.2655 0.0332 -0.0452 0.0356 -0.1146
Bias/constant: 0.3620

Evaluating model

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

MSE train = 0.0026
MSE test = 0.0020

R2 train = 0.9267
R2 test = 0.9300

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

Predicted y = 0.5329

End demo

The key calling code is:

LinearRegressor model = new LinearRegressor();
int nEpochs = model.Train(trainX, trainY);
Console.WriteLine("Done ");
Console.WriteLine("Used " + nEpochs + " epochs ");

The point is, the model is created and trained without any tuning. For use in a VIF system, the calling code can be simpler by ignoring the return value from Train():

LinearRegressor model = new LinearRegressor();
model.Train(trainX, trainY);

Good fun.

Implementing a machine learning regression system that just works without any manual input or tuning is difficult because there are so many things that can go wrong and can’t really be anticipated.

Right: I have never taken a selfie. This woman probably has second thoughts about taking selfies in the future after a more-then-close encounter with the camel.

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;

// adaptive LR: lr(t) = 0.01 / (t^0.25)
// stop: max_change_in_wts / max_weight "lt" 0.001
// max_iter = 10,000

namespace LinearRegressionSGDAdaptiveLR
{
  internal class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression" +
        " SGD training with adaptive LR and auto-exit ");

      // 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
      Console.WriteLine("\nCreating and training" +
        " Linear Regression model ");
      LinearRegressor model =
        new LinearRegressor();
      int nEpochs = model.Train(trainX, trainY);
      Console.WriteLine("Done ");
      Console.WriteLine("Used " + nEpochs + " epochs ");

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

      // 3. evaluate model
      Console.WriteLine("\nEvaluating model ");

      double accTrain = model.Accuracy(trainX, trainY, 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"));

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

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

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

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

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

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

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

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

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

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

    public int Train(double[][] trainX, double[] trainY,
      int maxEpochs = 10000, double exitTol = 0.001,
      double initRate = 0.01)
    {
      int n = trainX.Length;  
      int dim = trainX[0].Length;
      this.weights = new double[dim];

      // initialize weights and bias small rnd values
      double low = -0.01; double hi = 0.01;
      for (int i = 0; i "lt" dim; ++i)
        this.weights[i] = (hi - low) *
          this.rnd.NextDouble() + low;
      this.bias = (hi - low) *
          this.rnd.NextDouble() + low;

      int[] indices = new int[n];  // of train data
      for (int i = 0; i "lt" n; ++i)
        indices[i] = i;

      double[] prevWeights = new double[dim];
      for (int j = 0; j "lt" dim; ++j)
        prevWeights[j] = this.weights[j]; // for auto-exit
      
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        Shuffle(indices, this.rnd);
        //  lr(t) = 0.01 / (t^0.25)
        double lrnRate = 
          initRate / Math.Pow((double)(epoch+1), 0.25); 
        for (int i = 0; i "lt" n; ++i) // each train item
        {
          int ii = indices[i];
          double[] x = trainX[ii];
          double predY = this.Predict(x);
          double actualY = trainY[ii];
          for (int j = 0; j "lt" dim; ++j) // each weight
            this.weights[j] -= lrnRate *
              (predY - actualY) * x[j];
          this.bias -= lrnRate * (predY - actualY);
        }
        //if (epoch % (int)(maxEpochs / 5) == 0) // progress
        //{
        //  double mse = this.MSE(trainX, trainY);
        //  string s = "";
        //  s += "epoch = " + epoch.ToString().PadLeft(5);
        //  s += "  MSE = " + mse.ToString("F4").PadLeft(8);
        //  Console.WriteLine(s);
        //}

        // check for early-exit after each epoch
        // if max_change_in_wts / max_weights "lt" exit-tol
        double[] weightDeltas = new double[dim];
        for (int j = 0; j "lt" dim; ++j)
          weightDeltas[j] = 
            Math.Abs(prevWeights[j] - this.weights[j]);
        double maxChange = 0.0;
        for (int j = 0; j "lt" dim; ++j)
          if (weightDeltas[j] "gt" maxChange)
            maxChange = weightDeltas[j];
        double maxWeight = Math.Abs(this.weights[0]);
        for (int j = 0; j "lt" dim; ++j)
          if (Math.Abss(this.weights[j]) "gt" maxWeight)
            maxWeight = Math.Abs(this.weights[j]);
        if (maxWeight != 0.0 &&
          (maxChange / maxWeight) "lt" exitTol)
        {
          // Console.WriteLine("Early exit at epoch " +
          // epoch);
          return epoch;
        }

        // early exit didn't happen
        for (int j = 0; j "lt" dim; ++j)
          prevWeights[j] = this.weights[j];

      } // epoch
      return maxEpochs;

    } // Train

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

    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 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 * actualY))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

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

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

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

    public double R2(double[][] dataX, double[] dataY)
    {
      int n = dataX.Length;
      double sum = 0.0;

      for (int i = 0; i "lt" n; ++i)
        sum += dataY[i];
      double meanY = sum / n;

      double ssRes = 0.0; // sum squared residuals
      double ssTot = 0.0; // sum squared total
      for (int i = 0; i "lt" n; ++i)
      {
        double predY = this.Predict(dataX[i]);
        ssRes +=
          (dataY[i] - predY) * (dataY[i] - predY);
        ssTot +=
          (dataY[i] - meanY) * (dataY[i] - meanY);
      }
      return 1.0 - (ssRes / ssTot);
    }

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

    private static void Shuffle(int[] indices, Random rnd)
    {
      int n = indices.Length;
      for (int i = 0; i "lt" n; ++i)
      {
        int ri = rnd.Next(i, n);
        int tmp = indices[i];
        indices[i] = indices[ri];
        indices[ri] = tmp;
      }
    }

  } // class LinearRegressor

} // ns

Training data:

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

Test data:

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