Support Vector Regression From Scratch Using Python With SGD (SSGD) Training

Bottom line: I put together a demo of (kernelized) support vector regression that uses stochastic gradient descent (SGD) training. It works fine but training is relatively slow. Actually, the technique I used is stochastic sub-gradient descent (SSGD) but it’s common to refer to it as SGD.

The goal of a machine learning regression problem is to predict a single numeric value. Common regression techniques are linear regression, nearest neighbors regression, quadratic regression, kernel ridge regression (and the closely-related Gaussian process regression), neural network regression, random forest regression, and gradient boost regression. Each technique has many variations, and each technique has pros and cons.


The from scratch version (left) gives the same results as the scikit version (right), but the models have different weights.

Support vector regression (SVR) used to be popular in the late 1990s, for reasons which kind of baffle me. Kernel ridge regression is closely related to SVR and kernel ridge regression is easier to implement, easier to train, easier to interpret, and almost always gives better results than SVR (because KRR is easier to train). But, there are scenarios where SVR is required — typically legacy systems.

In almost all code libraries, SVR is trained using a form of quadratic programming or a strange algorithm called sequential minimal optimization (SMO). Both techniques are a nightmare to implement from scratch. In fact, SVR is so difficult to implement, to the best of my knowledge, every SVR library module I’ve seen relies on (is a wrapper around) a single C++ implementation called libsvm.

Kernelized SVR uses a kernel function, usually RBF (radial basis function). RBF requires a parameter usually called gamma (there’s a sigma version too). It’s possible to virtually reduce the number of training data items by driving their weights to zero — these are called the support vectors. In theory this leads to faster predictions, but in practice there is no increase in performance, except in rare scenarios.

Output of a demo of my from-scratch SVR:

Begin scratch SVR using SGD training

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

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

First three train y:
0.4840
0.1568
0.8054

Creating scratch Python SVR model
Setting gamma = 0.3000
Setting C = 0.999995
Setting epsilon = 0.003000
Setting lrn_rate = 0.0010
Setting max_epochs = 10000

Training SVR model using SGD
epoch =    0  |  MSE = 0.2874
epoch = 2000  |  MSE = 0.0000
epoch = 4000  |  MSE = 0.0000
epoch = 6000  |  MSE = 0.0000
epoch = 8000  |  MSE = 0.0000
Done

Model weights:
[-0.9999 -0.9999  0.0213 -0.6919  0.4950 . . .
  0.8400  0.6851  0.0000  0.0881  0.3400 . . .
 . . . 
 -0.9189  0.9999  0.0000 -0.9999 -1.0000 . . .
 -0.4469 -0.0001  0.0000  1.0000]
Number support vectors = 185

Train accuracy (0.10) = 0.9850
Test accuracy (0.10) = 0.9500

Train MSE = 0.0000
Test MSE = 0.0002

End demo

The output of a demo run using the scikit SVR module on the same data gives essentially the same results. I set the value of epsilon in my from-scratch implementation to 0.003 only to get identical results as the scikit version. The biggest practical downside to SVR is that it is very difficult to tune the RBF gamma, epsilon, and C parameters. My from scratch implementation adds learn_rate and max_epochs parameters to deal with too.

The scikit output:

Begin SVR using scikit

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 scikit SVR model
Setting gamma = 0.3000
Setting C = 0.999950
Setting epsilon = 0.0010
Done

Training scikit SVR model
Done

Mpdel weights:
[[-1.0000 -1.0000  0.1666 -0.9206  0.2359 -1.0000 . . .
   0.7708  1.0000  0.4633  0.0434 -0.0722 -0.3519 . . .
   . . . 
  -1.0000  1.0000  0.0383 -1.0000 -1.0000 -0.4875 . . .
  -0.3638 -0.5812  0.7083]]
Number model support vectors: [185]

Train accuracy (0.10) = 0.9850
Test accuracy (0.10) = 0.9500

Train MSE = 0.0000
Test MSE = 0.0002

End demo

I used a set of synthetic data that was generated by a neural network with random weights and biases. Each item has five predictor values. There are 200 training items and 40 test items.

The gamma parameter controls the behavior of the RBF function. The C (“complexity”) parameter is used for regularization to limit the magnitude of the model weights (there is one weight for each training/support item). The epsilon parameter defines how close a prediction must be to its target, in order to be ignored during SVR training.

I ran the scikit model first, to determine how many support vectors are generated for that model.

My from-scratch version does not use a bias term, which is usually OK as long as the data isn’t wildly skewed in some way. You can always normalize or center the training data if necessary (but it’s an annoying task). One significant downside to the SVR trained using SGD/SSGD idea is performance — it is much slower compared than the scikit version.

I’m not entirely satisfied with this implementation. When I get some time, I’ll refactor my code to add an explicit bias term to match the design of the scikit version.

A fascinating exploration.



One way to think about machine learning regression is that it’s a search for hidden patterns in data. And more abstractly, all of science is a sort of a search for hidden truth.

The covers of every issue of Playboy Magazine (except for the first issue in December 1953) has a company bunny logo somewhere. In many cases, the bunny logo is prominent and clearly visible. But some covers have the logo cleverly hidden.

Left: On the cover of the June 1991 issue, the logo is disguised as part of the straw thatching on the edge of the model’s hat. The logo is to the left of the ‘B’ in the “By James Jones” text.

Right: On the cover of the February 1994 issue, the logo is disguised as a reflection in the nail polish on the model’s right thumb.


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

# svr_sgd.py
# kernel support vector regression with SGD training

import numpy as np

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

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

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

def accuracy(model, data_X, data_y, pct_close):
  if data_X.size == 0: return 0.0
  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):
  if data_X.size == 0: return -1.0
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    pred_y = model.predict(x)[0]
    diff = pred_y - y
    sum += diff * diff
  return sum /n

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

class MySVR:
  def __init__(self, gamma=0.1, epsilon=0.1, C=1.0,
    lr=0.01, max_epochs=1000, seed=0):
    self.gamma = gamma
    self.epsilon = epsilon
    self.C = C
    self.lr = lr
    self.max_epochs = max_epochs
    self.weights = None
    self.train_X = None
    self.train_y = None
    self.rnd = np.random.RandomState(seed)

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

  def rbf(self, x1, x2):
    sum = 0.0
    for i in range(len(x1)):
      sum += (x1[i] - x2[i]) * (x1[i] - x2[i])
    result = np.exp(-1 * self.gamma * sum)
    return result

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

  def make_K(self, X):
    n = len(X)
    K = np.zeros((n,n))
    for i in range(0,n):
      for j in range(i,n):
        z = self.rbf(X[i], X[j])
        K[i,j] = z; K[j,i] = z
    return K

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

  def fit(self, X, y):
    self.train_X = X.copy()
    self.train_y = y.copy()
    n, dim = self.train_X.shape

    # init weights
    self.weights = np.zeros(n)
    lo = -0.10; hi = 0.10
    for i in range(n):
      self.weights[i] = (hi - lo) * self.rnd.random() + lo

    K = self.make_K(self.train_X)  # lookup fast predicts
    lamda = 1.0 / self.C
    freq = self.max_epochs // 5
    
    for epoch in range(self.max_epochs):
      indices = self.rnd.permutation(n)
      for i in indices:
        y_pred = np.dot(K[i], self.weights)  # fast
        # y_pred = self.predict_one(X[i])  # slow!
        error = y_pred - self.train_y[i]

        grad_reg = lamda * self.weights[i]
        if error "gt" self.epsilon:
          grad_loss = 1.0
        elif error "lt" -self.epsilon:
          grad_loss = -1.0
        else:
          grad_loss = 0.0  # ignore inside epsilon tube
        self.weights[i] -= self.lr * (grad_reg + grad_loss)

      if epoch % freq == 0:
        m = mse(self, self.train_X, self.train_y)
        print("epoch = %4d  |  MSE = %0.4f " % (epoch,m))

    return  # all done

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

  def predict_one(self, x):
    # helper for predict(X)
    sum = 0.0
    for i in range(len(self.weights)):
      sum += self.weights[i] * self.rbf(x, self.train_X[i])
    return sum

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

  def predict(self, X):
    # X is a matrix of input vectors (scikit API)
    preds = []
    for i in range(len(X)):
      py = self.predict_one(X[i])
      preds.append(py)
    return np.array(preds)

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

  def get_supp_idxs(self):
    result = []
    for i in range(len(self.weights)):
      # a nearly-zero wt is associated with a supp vector
      if np.abs(self.weights[i]) "gt" 1.0e-5:
        result.append(i)
    return result

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

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

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

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

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

  # Creating scikit SVR model
  # Setting gamma = 0.3000
  # Setting C = 0.999950
  # Setting epsilon = 0.0010
  # Number model support vectors: [185]
  # Train accuracy (0.10) = 0.9850
  # Test accuracy (0.10) = 0.9500  
  # Train MSE = 0.0000
  # Test MSE = 0.0002

  # create and train model
  print("\nCreating scratch Python SVR model ")
  gamma = 0.30

  # smaller epsilon == fewer ignored == more supp vecs
  # larger epsilon == more ignored == fewer supp vecs
  epsilon = 0.003

  C = 0.999995
  lr = 0.001
  max_epochs = 10000

  print("Setting gamma = %0.4f " % gamma)
  print("Setting C = %0.6f " % C)
  print("Setting epsilon = %0.6f " % epsilon)
  print("Setting lrn_rate = %0.4f " % lr)
  print("Setting max_epochs = " + str(max_epochs))

  print("\nTraining SVR model using SGD ")
  model = MySVR(gamma=gamma, epsilon=epsilon,
    C=C, lr=lr, max_epochs=max_epochs, seed=1)
  model.fit(train_X, train_y)
  print("Done ")

  print("\nModel weights: ")
  print(model.weights)

  supp_vec_idxs = model.get_supp_idxs()
  print("Number support vectors = " + \
    str(len(model.get_supp_idxs())))
    
  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)

  print("\nEnd demo ")

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

if __name__ == "__main__":
  main()

The scikit SVR demo program.

# svr_scikit.py
# scikit-learn SVR module

import numpy as np
from sklearn.svm import SVR

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

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

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

def accuracy(model, data_X, data_y, pct_close):
  if data_X.size == 0: return 0.0
  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):
  if data_X.size == 0: return -1.0
  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 main():
  print("\nBegin SVR using scikit ")

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

  # create and train model
  # SVR(*, kernel='rbf', degree=3, gamma='scale',
  # coef0=0.0, tol=0.001, C=1.0, epsilon=0.1, 
  # shrinking=True, cache_size=200, verbose=False,
  # max_iter=-1)

  print("\nCreating scikit SVR model ")
  gamma = 0.30
  epsilon = 0.001
  C = 0.99995

  print("Setting gamma = %0.4f " % gamma)
  print("Setting C = %0.6f " % C)
  print("Setting epsilon = %0.4f " % epsilon)
  model = SVR(kernel='rbf', gamma=gamma, C=C, 
    epsilon=epsilon)
  print("Done ")

  print("\nTraining scikit SVR model ")
  model.fit(train_X, train_y)
  print("Done ")

  print("\nMpdel weights: ")
  print(model.dual_coef_)

  print("Number model support vectors: " + \
    str(model.n_support_))

  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)

  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
Posted in Machine Learning, Scikit | Leave a comment

Machine Learning Myths: The Bias Term in Kernel Ridge Regression Models and in Support Vector Regression Models

In spite of decades of research in machine learning, there are still dozens of pieces of common, but just plain incorrect, pieces of knowledge floating about on the Internet. I came across one recently. Does a kernel ridge regression model need a bias term or not?

Let me cut to the chase and say that the answer is that kernel ridge regression models do not need a bias term. But many Internet and AI resources claim that a bias term is absolutely needed. For example, an AI told me today:

“The bias term (intercept) in a kernel ridge regression model acts as an offset to shift the regression line or hyperplane up or down. Without it, the model is geometrically constrained to pass exactly through the origin (f(x) = 0 when (x = 0), which would drastically hurt predictive accuracy.”

Wrong.

The correct logic is:

“Kernel ridge regression works well without an explicit bias term because the kernel trick implicitly maps data into a high-dimensional space where “non-centered” data is inherently handled by the flexibility of the kernel, and the similarity computations act as a localized center.

Many kernels (including the radial basis function) measure the localized similarity between points. The prediction is essentially computed as a weighted sum of similarities to the training vectors. Even if raw training data is not centered, this process simply shifts the “humps” of the kernel evaluations around, allowing the model to naturally fit the true offsets.”

The scikit KernelRidge module does not introduce a bias term into its models.



Output of a demo:

Begin demo

Generating 40 rows dummy data
Done

First three X:
[[1.0976]
 [1.4304]
 [1.2055]]

First three y:
[0.8013 0.7921 0.8992]

Creating and training SVR model
Done

Model bias = 0.6412
Model R2 score = 0.5289

Creating and training KRR model
Done

Model has no bias
Model R2 score = 0.6373

End demo

The moral of the story is that it’s prudent to be wary of AI-generated information. AI always presents information with a sense of complete confidence. Humans are conditioned to believe just about anything, from anyone or anything, when it is presented with an air of absolute confidence. (Including information in blog posts like this one).



Machine learning myths have relatively low impact. But some myths have very big consequences. One common myth that I hear, primarily from young people in their 20s and 30s, is that intelligence has no genetic component, and therefore is not inherited.

There is overwhelming scientific evidence that intelligence is at least 50% inherited, and most likely closer to about 80% inherited. Put simply, smart parents produce smart children, and low-IQ parents produce low-IQ children.

Put more simply, “Stupid breeds stupid.”

I don’t understand why intelligence-inheritability is aggressively censored by virtually all of mainstream media, and is almost 100% censored in academia.

When I was a college student at UC Irvine, I worked at Disneyland in the evenings and the weekends. Everyone knows that the Park’s motto is, “The Happiest Place on Earth”. I’m going to guess that these brawling families at Disneyland are not examples of inherited intelligence at its finest, and that they are not having an especially Happy Time. And they probably had even less fun after they were all arrested.


Demo program:

# svr_bias_scikit.py

import numpy as np
from sklearn.svm import SVR
from sklearn.kernel_ridge import KernelRidge

import random
import numpy as np
random.seed(0)
np.random.seed(0)

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

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

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

print("\nBegin demo ")

print("\nGenerating 40 rows dummy data ")
X = 2 * np.random.rand(40, 1)
y = np.sin(X).ravel() + \
  np.random.normal(0, 0.1, X.shape[0])
print("Done ")

print("\nFirst three X: ")
print(X[0:3])

print("\nFirst three y: ")
print(y[0:3])

gamma = 0.02

print("\nCreating and training SVR model ")
svr = SVR(gamma=gamma)
svr.fit(X, y)
print("Done ")

if hasattr(svr, "intercept_"):
  svr_bias = svr.intercept_[0]
  print("\nModel bias = %0.4f " % svr_bias)
else:
  print("\nModel has no bias ")

svr_r2 = svr.score(X, y)
print("Model R2 score = %0.4f " % svr_r2)

print("\nCreating and training KRR model ")
krr = KernelRidge(gamma=gamma)
krr.fit(X, y)
print("Done ")
if hasattr(krr, "intercept_"):
  krr_bias = krr.intercept_[0]
  print("\nModel bias = %0.4f " % krr_bias)
else:
  print("\nModel has no bias ")

krr_r2 = krr.score(X, y)
print("Model R2 score = %0.4f " % krr_r2)

print("\nEnd demo ")
Posted in Machine Learning, Scikit | Leave a comment

Linear Regression with SGD Adaptive Learning Rate and Auto-Exit Using C#

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.

On the other hand, in life, there are some scenarios where the “What could possible go wrong” question isn’t too difficult to answer.

Left: The California Alligator Farm was a major attraction in Los Angeles from 1907 to 1953, and then after a move to Buena Park in 1953, was a very popular spot until the farm closed in 1984. I grew up in Anaheim and Fullerton, two cities adjacent to Buena Park. My family went to the Alligator Farm often in the 1960s.

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
Posted in Machine Learning | Leave a comment

My Top Ten Favorite Science Fiction Movies of All Time

I am a big fan of science fiction movies. Here’s a list of my top ten favorites (plus one extra) of all time. My main criterion is, if I were going to a remote Arctic research station and could take only 10 science fiction movies, which 10 would I bring? Every movie on this list has a great story, great special effects, and good-or-excellent acting. Listed by year released.


1. Godzilla (1954/1956) – I like both the original 1954 Japanese version and the 1956 American adaptation. Unlike later Godzilla movies, the original is deadly serious. The special effects are amazing. The scene where Godzilla first appears over a hill on Odo Island, terrified me as a young man and gave me nightmares for many years afterwards. My grade = A+.


2. Forbidden Planet (1956) – In the 23rd century, the crew of the C57-D travel to planet Altair IV to find out what happened to an expedition 20 years earlier. The only survivors are Dr. Morbius and his 20-years-old daughter Altaira. The crew is menaced by a seemingly unstoppable monster made of energy. Great story, great acting, great special effects, ground-breaking sound effects. I’m always baffled why this movie is usually left off most other top-ten sci-fi movie lists. My grade = A.


3. Star Wars (1977) – It’s almost impossible to overstate what a profound impact this movie had when it was released. It’s arguably the most famous and influential science fiction movie in history. I think there are better movies than “Star Wars on this list, but “Star Wars” is absolutely a top ten entry. My grade = A.


4. Alien (1979) – The crew of the tramp spacecraft Nostromo investigates a crashed alien spaceship on a creepy planet. Bad idea. The scene where the alien pops out of a crewman’s chest is one of the most famous in science fiction movie history. I don’t like most sci-fi horror films, but this movie is the one exception. I vividly remember watching this movie on the day it was released in 1979, along with some of my friends who worked at Disneyland in Anaheim with me. An absolutely terrifying film. My grade = A.


5. Jurassic Park (1993) – Everyone I know has reached a Jurassic movie fatigue with all the Jurassic sequels, but the first movie in the series is great and is the by far best. The special effects were absolutely astonishing in 1993 — an incredible advance over anything previously seen. My grade = A.


6. Starship Troopers (1997) – Many professional movie critics don’t like this movie, but it easily makes my top ten of all time. Sometime in the not-too-distant future, Earth wages war against a very nasty alien insectoid species. Things go bad for Earth early on, but humans prevail in the end. Brilliant special effects and exciting action. My grade = A+.


7. The Fifth Element (1997) – Another movie that many professional critics don’t like very much. But I love this movie. Wildly creative story (even if it doesn’t make a whole lot of sense) and wildly creative details. I can understand why some people don’t like “The Fifth Element”, but I’m a big fan. My grade = A+.


8. The Matrix (1999) – This movie makes almost all top-ten sci-fi movie lists and I agree. Thomas Anderson, aka Neo, discovers that robots and AI have enslaved all humans, putting people in pods and making them bio-batteries. This movie has so many nuances, I enjoyed watching it several times to pick up details I hadn’t noticed before. My grade = A+.


9. Inception (2010) – A team of high-tech agents steals industrial information by entering the dreams and subconsciousnesses of their targets. Dreams, within dreams, within dreams. A very imaginative movie with great special effects. My grade = A.


10. Dune: Part One (2021) – An epic story where the royal family House of Atreides is menaced by the House of Harkonnen, the evil Harkonnen Saraukar troops, and the backing of the Emperor. Based closely on one of my all-time favorite science fiction novels. I liked the 1984 version of “Dune” but the 2021 version is a rare example of a sci-fi movie remake that is significantly better than the original. My grade = A.


11. Project Hail Mary (2026) – In 2032, scientist-teacher Ryland Grace must travel to star Tau Ceti to find a way to stop mysterious astrophages from destroying Earth’s sun. He meets an alien who is on the same mission to save his planet Eridia. Together they succeed. Great story, great acting, great special effects. A new movie on this list. It is tentatively joining while I wait to see if it stands the test of time. My grade = A+.



12. Dark City (1998) – A man wakes up in a strange city and has no memory of who he is. There is only constant night, and mysterious, menacing strangers. It turns out the man is in an alien experiment in space. A clever-twist, happy, ending. This movie was on my top ten list for years, but was bumped off the list by “Inception”. My grade = A.


Honorable Mentions – Movies that Barely Miss My Top Ten


“The Thing from Another World “(1951) – An unfriendly alien crash lands near an Arctic research station. My grade = A-.

“Invaders from Mars” (1953) – A Martian advance party to prepare a full-scale invasion. Has top 10 sound effects. The path to the sandpit terrified every person I know who saw this movie as a child (including me). My grade = A-.

“War of the Worlds” (1953) – A full-scale Martian invasion. Academy Award for Special Effects. My grade = B+.

“Gog” (1954) – Strange deaths at a secret underground laboratory that has a robot named Gog — that has an attached flamethrower, because all research robots need a flamethrower. My grade = B+.

“Quatermass 2” aka “Enemy from Space” (1957) – Aliens use parasites to control UK villagers to force them to build a plant to make food and an artificial environment to prepare for an invasion. Has top 10 sound effects. My grade = B+.

“Bladerunner” (1982) – In a futuristic 2019, an agent (who may be an android himself) tracks down rogue androids. My grade = B+.

“Predator Badlands” (2025) – A young alien must prove his manhood (well, I suppose it’s “alien-hood”) by capturing the most dangerous creature known on the Death Planet. My grade = A-.



Dishonorable Mentions: Movies On Most Top-Ten Lists But Ones I Just Don’t Like At All


“The Planet of the Apes” (1968) – Ridiculous and annoying. My grade = D.

“Any of the Zillions of Other Ape Movies” – All of them terrible. Apes belong in zoos or politics, not in sci-fi movies. My grade = D.

“2001: A Space Odyssey” (1968) – Mind-boggling nonsense that appeals only to pseudo-intellectuals who are too insecure to admit they don’t understand it, and so they won’t say they don’t like it. My grade = D+.

“Close Encounters of the Third Kind” (1977) – A close encounter with boredom. My grade = D.

“E.T. the Extra-Terrestrial” (1982) – Annoying kids, annoying alien, annoying movie. My grade = D-.

“Avatar” (2009) – I should like it, but something about this movie just didn’t work for me. My grade = C-.


Posted in Top Ten | Leave a comment

“Support Vector Regression with SGD Training Using C#” in Visual Studio Magazine

I wrote an article titled “Support Vector Regression with SGD Training Using C#” in the July 2026 issue of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2026/07/01/support-vector-regression-with-sgd-training-using-csharp.aspx.

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

There are approximately a dozen common regression techniques. Each technique has pros and cons. A technique that sometimes produces highly accurate predictions for specific types of data is called support vector regression (SVR).

The term “support vector regression” by itself is ambiguous, because there are two types of SVR: linear SVR and kernelized SVR. Linear SVR is rarely used so the term “SVR” by itself usually means kernel SVR.

My article presents a demo of kernel SVR, trained using a variation of stochastic gradient descent (SGD) called stochastic sub-gradient descent (SSGD). Note: Even though SSGD is technically a variation of SGD, the two techniques are so similar that SSGD is often called SGD.

Support vector regression (SVR) predicts using a kernel function, usually RBF, that computes similarity between two data items. During training, some of the training items are determined to be irrelevant to some extent, and they are removed. The training items that are left are called the support vectors.

The key parts of the demo output are:

Creating SVR object
Setting RBF gamma = 0.3000
Setting epsilon = 0.007500
Setting C = 1.00
Setting lrnRate = 0.0010
Setting maxEpochs = 5000
Setting tol = 0.000100

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

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

Model bias = 0.4030
Number supp vectors = 194

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

The first versions of support vector regression in the 1980s required quadratic programming optimization training, which is very complicated, slow, and doesn’t scale well to large datasets. The sequential minimal optimization (SMO) training algorithm was developed in the late 1990s. SMO is fast but the algorithm is complicated, and very difficult to correctly implement. The SGD (actually stochastic sub-gradient descent, SSGD) training technique presented in the article is by far the simplest training algorithm, and often works well in practice because processing one data item at a time introduces a form of implicit regularization, which helps the SVR model to predict new, previously unseen data.

Support vector regression had a brief surge of popularity in the late 1990s and early 2000s. However, data scientists realized that the closely related kernel ridge regression (KRR) has several significant advantages over SVR, and so the use of SVR declined to the point where it is not used very much today. SVR is more difficult to implement than KRR, SVR is much more difficult to tune than KRR (KRR can use true SGD, which is easier to tune than SVR sub-gradient descent), and SVR often gives slightly worse prediction accuracy than KRR (due mostly to the difficulty in parameter tuning). That said, there are some problem domains, such as biology and chemistry, where kernel SVR is often used and is highly effective.



I don’t collect things. I collect weird experiences. A few years ago, I noticed that when I have my TV closed captioning turned on, a surprising number of animals “chitter”. Now I’m always on the alert for “chittering”.

Racoons are notorious chitterers.

Left: “101 Dalmatians” (1996) – You know the story. The evil Cruella de Vil wants to make a coat out of the pelts of Dalmatian puppies. The 1961 animated version is excellent. This live-action version is mediocre, maybe a C+ grade from me.

Center: “Pocahontas” (1995) – Disney had several great animated films in the 1990s — this wasn’t one of them. Beautiful animation but weak, woke story. Grade = C.

Right: “Holland” (2025) – I had hopes for this movie because it stars Nicole Kidman who is a terrific actress. Nope. A woman discovers she married a serial murderer. Huh? A racoon in an alley trash can was the highlight of the movie. Grade = C-.


Posted in Machine Learning | Leave a comment

Support Vector Regression Trained With Particle Swarm Optimization: Failure

I rarely post failed attempts, but a recent failure illustrates the point that most failures provide a lesson, and that lesson can lead to a success. Briefly, I implemented support vector regression (SVR) training using particle swarm optimization (PSO). The experiment was a failure in the sense that the trained SVR model predicted at only 35% accuracy, while an SVR model trained using standard quadratic programming techniques predicts at 95% accuracy.

Training an SVR model is a good challenge because the underlying optimization problem (epsilon-insensitive loss with L2 regularization) is not Calculus-differentiable, and so stochastic gradient descent (SGD) cannot be used (but a variant called stochastic sub-gradient descent is possible). Additionally, SVR training has no matrix-inverse closed form solution. Therefore, one of the few practical ways to train an SVR model is to use quadratic programming — a hideous technique.



A screenshot of one of the dozens of failed attempts.


These facts are the motivation for looking at training using PSO. In theory, PSO can solve any optimization problem, including training an SVR model. But in practice, PSO often just doesn’t work.

Here’s the output of the most successful of my many failed attempts to train an SVR model using PSO:

Begin Support Vector Regression with particle
 swarm optimization training

Loading train (200) and test (40) from file
Done

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

First three train y:
  0.4840
  0.1568
  0.8054

Setting kernel SVR parameters:
C = 10.00
epsilon = 0.01
gamma = 0.10

Setting particle swarm training parameters:
numParticles = 100
maxIter = 1000
lrnRate = 0.0500

Starting PSO training
Begin prelim pseudo-SGD pass
Finished
iteration =        0  ep-insensitive loss = 116.4313
iteration =      100  ep-insensitive loss =  69.2760
iteration =      200  ep-insensitive loss =  64.7903
iteration =      300  ep-insensitive loss =  64.0210
iteration =      400  ep-insensitive loss =  63.3281
iteration =      500  ep-insensitive loss =  62.5250
iteration =      600  ep-insensitive loss =  62.0185
iteration =      700  ep-insensitive loss =  60.6460
iteration =      800  ep-insensitive loss =  59.9317
iteration =      900  ep-insensitive loss =  58.6440
Number primal support vectors = 161
Done

Model wts:
  -0.0002  -0.0822   0.1847   0.0870   0.1146   0.0000
   0.1390   0.0038   0.0000   0.0000   0.1325  -0.0112
   0.0389  -0.2398   0.0000  -0.0769  -0.1857   0.0000
   0.1819  -0.2268  -0.0253  -0.0337   0.0000   0.1092
   0.0462  -0.0352   0.0000  -0.0509  -0.0779   0.0000
   0.0000   0.0390  -0.0398   0.0000  -0.0097   0.0582
  -0.1024   0.0033  -0.0941   0.0000  -0.0024  -0.0111
  -0.0804  -0.1082   0.0000  -0.0585  -0.0311  -0.0228
  -0.0530  -0.0708  -0.1301   0.0085  -0.0487  -0.1630
   0.0159   0.0807  -0.2094  -0.0207   0.0701   0.1473
   0.0000   0.0000   0.0869   0.0000   0.0000  -0.0769
   0.0000   0.4215   0.0000  -0.1014  -0.1130   0.0000
   0.0430  -0.0371  -0.0423  -0.0721  -0.0289   0.2837
  -0.0245   0.0668   0.0000  -0.0430  -0.2214   0.2056
   0.0712   0.0526   0.0297   0.0198   0.0850  -0.1091
  -0.1088  -0.0209   0.0678  -0.1045  -0.1739   0.0764
   0.1685   0.0059   0.1136   0.0000   0.0078   0.0000
   0.0000   0.0948   0.0495   0.0000   0.0147   0.1114
  -0.0101   0.2158  -0.1645  -0.0221  -0.1714  -0.0395
  -0.1103   0.0000   0.1159   0.0264  -0.0595  -0.0855
   0.0000  -0.0742  -0.0558  -0.0898   0.0239  -0.0950
  -0.0596   0.0355  -0.0018   0.0984  -0.0520   0.0229
  -0.1841  -0.0647   0.1025   0.0611  -0.0307  -0.0287
   0.0000  -0.0200   0.0000  -0.1502  -0.0338   0.0366
   0.0673   0.0665   0.0000   0.0304  -0.0521   0.0727
   0.3018  -0.0553   0.0000   0.0776   0.0000   0.0419 
   0.0420   0.0000   0.0219   0.0758   0.0000   0.0000
   0.0000   0.0361   0.0955  -0.0149   0.0182  -0.0562
   0.0985  -0.0833   0.0774  -0.1483  -0.0067   0.0000
  -0.0732  -0.2078  -0.0735  -0.0280  -0.0254   0.0312
  -0.1959  -0.0398   0.0000  -0.0135   0.0000  -0.0047
  -0.0112   0.0330   0.0121  -0.0104   0.0837   0.2852
   0.0394   0.4294   0.0000   0.2324   0.0263   0.1959
   0.1695   0.0153

Computing model accuracy

Train acc (within 0.10) = 0.3500
Test acc (within 0.10) = 0.3500

End demo

In SVR, each training item has an associated weight that is used when making a prediction. The demo training data has 200 items so there are 200 weights.

Weights that are zero have no effect. The non-zero weights are associated with training items called support vectors. In a non-demo scenario, I would remove the training items associated with zero weights, and remove the zero-weights.

In previous failed attempts, the PSO never got started. Therefore for the attempt shown, I ran the data through a preliminary pseudo stochastic gradient descent phase to get a good starting point for the PSO. This worked, but then a new problem arose. The PSO training converges quickly to a solution with loss about 65.00 and then stalls out.

This is an example of the exploration-exploitation issue in machine learning. This demo has too much exploitation (looking at new possible solutions that are very close to the current best known solution), and not enough exploration (looking at completely new solutions, to avoid getting stuck).

Therefore, the lesson learned is for me to discard the initial pseudo-SGD training, and modify the PSO training to increase the exploration part of the code.



Software design failure is one thing. Hotel name failure is another. Here are two foreign hotel signs that don’t strike a good note in American English.


Posted in Machine Learning | Leave a comment

Linear Ridge Regression From Scratch Using C# to Sync With the scikit Ridge Module

Linear ridge regression (often called just ridge regression, but not to be confused with kernel ridge regression) uses mathematical L2 regularization to prevent model weights from becoming too large, which in turn limits model overfitting. Whew! There’s a lot of ideas in that sentence.

Behind the scenes, L2 regularization penalizes the sum of the squared model weights, by conceptually adding that sum to the error function, which in turn modifies the gradient, which in turn modifies the weight updates. The ideas are relatively simple, but implementation is surprisingly subtle and tricky.

Ridge regression requires an alpha constant that controls the amount of L2 penalty. There are many ways to implement L2 regularization and all of them influence the meaning of the alpha constant.

A long time ago, I examined the Ridge module in the Python language scikit-learn library. The documentation specifies that alpha is “Constant that multiplies the L2 term, controlling regularization strength. Alpha must be a non-negative float. When alpha = 0, the objective is equivalent to ordinary least squares.”

This confused me because this definition did not correspond to any of the L2 implementations I had ever seen or implemented. So, I set out to understand how the scikit Ridge module implemented L2. After a lot of work I determined that scikit normalizes the L2 error term by dividing by the number of training items. In code:

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

  double err = predY - actualY;
  for (int j = 0; j "lt" dim; ++j) // each weight
  {
    // normed by train size n to sync with scikit
    double dw = err * x[j] + 
      (alpha * this.weights[j] / n); 
    this.weights[j] -= lrnRate * dw;
  }
  this.bias -= lrnRate * err; // no L2 on bias
}

The divide by n term was the key. Here’s the output of one run of my C# demo:

Begin ridge regression with SGD training using C#

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 ridge regression model
Done

Training model using SGD
Setting lrnRate = 0.0010
Setting maxEpochs = 1000
Setting L2 alpha = 0.5000
epoch =     0  MSE =   0.1132
epoch =   200  MSE =   0.0026
epoch =   400  MSE =   0.0026
epoch =   600  MSE =   0.0026
epoch =   800  MSE =   0.0026
Done

Weights/coefficients:
-0.2637  0.0331  -0.0453  0.0355  -0.1139
Bias/constant: 0.3619

Evaluating model

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

MSE train = 0.0026
MSE test = 0.0020

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

End demo

These C# results were identical (subject to rounding errors) to the output of a scikit Ridge demo.

My C# implementation uses standard stochastic gradient descent (SGD), while the scikit Ridge module uses an annoying complicated variation of SGD called SAG (stochastic average descent). Additionally, the Ridge module uses automatic stopping (when model weights don’t significantly change), and a programmatically generated learning rate that is based on the maximum of the sum of the squared values in each row of training data.

Note that for simplicity, I did not scale the training data to mean = 0 and variance = 1, which is more or less required in a non-demo scenario, so that all model weights are penalized equally.

A very interesting and satisfying exploration.



Optimization algorithms such as SGD and SAG are key components of a machine learning regression system.

The aircraft used by an airline is a key component. The Lockheed Constellation was the last major propellor-engine passenger plane before the jet aircraft Douglas DC-8 and Boeing 707 took over in 1959, almost overnight. There were about 850 Constellations produced from 1943-1958. The point is that new technologies can render old technologies obsolete almost instantly.


Demo C# program. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols.

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

// ridge linear regression with SGD training
// replicate scikit Ridge results

namespace LinearRegressionRidge
{
  internal class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin ridge regression" +
        " with SGD training using C# ");

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

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

      string testFile =
        "..\\..\\..\\Data\\synthetic_test_40.txt";
      double[][] testX = MatLoad(testFile, colsX, ',', "#");
      double[] testY = MatToVec(MatLoad(testFile,
        new int[] { colY }, ',', "#"));
      Console.WriteLine("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 model
      Console.WriteLine("\nCreating ridge regression model ");
      LinearRegressor model = new LinearRegressor();
      Console.WriteLine("Done");

      // 3. train model using SGD
      Console.WriteLine("\nTraining model using SGD ");
      double lrnRate = 0.001;
      int maxEpochs = 1000;
      double alpha = 0.5;
      Console.WriteLine("Setting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpochs = " +
        maxEpochs);
      Console.WriteLine("Setting L2 alpha = " +
        alpha.ToString("F4"));
      model.TrainSGD(trainX, trainY, lrnRate,
        maxEpochs, alpha);
      Console.WriteLine("Done ");

      // 4. examine 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"));

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

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

      // 7. TODO: save model weights, bias to text file

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

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

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

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

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

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

  public class LinearRegressor
  {
    public double[] weights;
    public double bias;
    private Random rnd; // for SGD training

    public LinearRegressor(int seed = 0)
    {
      this.weights = new double[0]; // quasi-null
      this.bias = 0.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 double Accuracy(double[][] dataX, double[] dataY,
      double pctClose)
    {
      int numCorrect = 0; int numWrong = 0;
      for (int i = 0; i "lt" dataX.Length; ++i)
      {
        double actualY = dataY[i];
        double predY = this.Predict(dataX[i]);
        if (Math.Abs(predY - actualY) "lt" 
          Math.Abs(pctClose * actualY))
          ++numCorrect;
        else
          ++numWrong;
      }
      return (numCorrect * 1.0) / (numWrong + numCorrect);
    }

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

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

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

    public void TrainSGD(double[][] trainX, double[] trainY,
      double lrnRate, int maxEpochs, double alpha)
    {
      int n = trainX.Length;
      int dim = trainX[0].Length;
      this.weights = new double[dim];
      double lo = -0.01; double hi = 0.01;
      for (int i = 0; i "lt" dim; ++i)
        this.weights[i] = (hi - lo) * 
          this.rnd.NextDouble() + lo;
      this.bias = (hi - lo) * this.rnd.NextDouble() + lo;

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

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

          double err = predY - actualY;
          for (int j = 0; j "lt" dim; ++j) // each weight
          {
            // normed by train size to sync with scikit
            double dw = err * x[j] + 
              (alpha * this.weights[j] / n); 
            this.weights[j] -= lrnRate * dw;
          }
          this.bias -= lrnRate * err; // no L2 on bias
        }
        if (epoch % (int)(maxEpochs / 5) == 0) // 5 times
        {
          double mse = this.MSE(trainX, trainY);
          string s1 = "epoch = " + 
            epoch.ToString().PadLeft(5);
          string s2 = "  MSE = " + 
            mse.ToString("F4").PadLeft(8);
          Console.WriteLine(s1 + s2);
        }
      }
    } // TrainSGD()

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

    private static void Shuffle(int[] indices, Random rnd)
    {
      // helper for TrainSGD()
      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

Demo scikit program.

# ridge_regression_scikit.py
# (linear) ridge regression demo

import numpy as np
from sklearn.linear_model import Ridge

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

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

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
  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):
    actual_y = data_y[i]
    pred_y = model.predict(data_X[i].reshape(1, -1))[0]
    diff = actual_y - pred_y
    sum += diff * diff
  return sum /n

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

print("\nBegin (linear) ridge regression using scikit ")

print("\nLoading 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("\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])

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

print("\nCreating ridge model with SAG training ")
alpha = 0.5
print("Setting L2 alpha = %0.4f " % alpha)
model = Ridge(alpha=alpha, solver='sag')
print("Done ")

print("\nTraining model with SGD ")
print("Using default training params ")
model.fit(train_X, train_y)
print("Done ")

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

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

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

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

print("\nEnd demo ")

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
Posted in Machine Learning, Scikit | Leave a comment

Linear Ridge Regression From Scratch With Stochastic Average Gradient (SAG) Training Using C#

One evening, I just couldn’t sleep, so I decided to implement linear ridge (aka L2 regularization) regression using the stochastic average gradient (SAG) training algorithm, from scratch, using the C# language. SAG is a variation of basic stochastic gradient descent (SGD).

Before I go any further, let me note that my initial thought was to implement standard linear regression using SAG, but my efforts were not very successful, and so I had to use ridge regression instead. It seems that SAG is designed to work with L2, but the details aren’t clear to me.

Bottom line: Although my experiments were not conclusive, SAG was more complicated than SGD to implement, SAG was harder than SGD to tune (learning rate and exit exit stop tolerance), and SAG uses more memory than SGD. I could not detect any speed/performance differences between SAG and SGD. In short, I could see no reason to use SAG instead of SGD (in my to-be-sure very limited experiments).

For my experiments, I used one of my standard synthetic datasets. The data looks like:

-0.1660,  0.4406, -0.9998, -0.3953, -0.7065,  0.4840
 0.0776, -0.1616,  0.3704, -0.5911,  0.7562,  0.1568
-0.9452,  0.3409, -0.1654,  0.1174, -0.7192,  0.8054
 0.9365, -0.3732,  0.3846,  0.7528,  0.7892,  0.1345
. . .

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

After displaying the first few lines in the training data, the output of the demo program is:

Begin C# linear regression SAG training (with early exit and L2)

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

Setting lrnRate = 0.0010
Setting maxEpohcs = 500
Setting stop tol = 0.000100
Setting alpha = 0.5000

Creating and training Linear Regression model using SAG
epoch =     0  MSE =   0.1150
epoch =   100  MSE =   0.0026
Early exit at epoch 131
Done

Weights/coefficients:
-0.2637 0.0332 -0.0454 0.0355 -0.1139
Bias/constant: 0.3618

Evaluating model

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

MSE train = 0.00258
MSE test = 0.00195

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

End demo

The model prediction accuracy is poor because the synthetic data has a complex non-linear structure.

I validated my demo against the Python language scikit-learn library Ridge module, which has an option to use SAG (along with more common algorithms such as L-BFGS and Cholesky). The scikit results were the same, to 4 decimal places, as my C# demo results.

I get the feeling that SAG was a research investigation in search of a problem, rather than an engineering project to solve a specific problem. During my years working at Microsoft Research, I saw this research-oriented strategy many times.



There are many different algorithms that can be used to train a linear regression model. Each algorithm has a different look and feel.

I learned to read in the 1950s and 1960s from comic books. Superman was one of my favorites. Most of the cover art in the 1950s was done by artist Wayne Boring (1905-1987). Most of the cover art in the 1960s was done by the famous Curt Swan (1920-1996). Artist Al Plastino (1921-2013) chipped quite a few covers in the 1950s, 60s, and 70s.

Left: “Superman” #111, February 1957. Cover art by Al Plastino.

Center: “Superman” #106, July 1956. Cover art by Wayne Boring.

Right: “Superman” #149, November 1961. Cover art by Curt Swan.


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

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

namespace LinearRegressionSAG
{
  internal class LinearRegressionSAGProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression" +
        " SAG training (with early exit and L2) ");

      // 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
      double lrnRate = 0.001;
      int maxEpochs = 500;
      double stopTol = 0.0001; // early training exit
      double alpha = 0.5; // L2 regularization, "ridge"
      //int seed = 0;
      Console.WriteLine("\nSetting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpohcs = " +
        maxEpochs);
      Console.WriteLine("Setting stop tol = " +
        stopTol.ToString("F6"));
      Console.WriteLine("Setting alpha = " +
        alpha.ToString("F4"));

      Console.WriteLine("\nCreating and training" +
        " Linear Regression model using SAG ");
      LinearRegressor model =
        new LinearRegressor();
      model.TrainSGD(trainX, trainY, lrnRate,
        maxEpochs, stopTol, alpha);
      Console.WriteLine("Done ");

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

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

      double accTrain = model.Accuracy(trainX, trainY, 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("F5"));
      double mseTest = model.MSE(testX, testY);
      Console.WriteLine("MSE test = " +
        mseTest.ToString("F5"));

      // 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 TrainSGD(double[][] trainX,
      double[] trainY, double lrnRate, int maxEpochs,
      double stopTol, double alpha)
    {
      int n = trainX.Length;  int dim = trainX[0].Length;
      this.weights = new double[dim];

      // initialize weights and bias
      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;

      // SAG memory
      double[][] gradMem = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        gradMem[i] = new double[dim];
      double[] totGrad = new double[dim];
      double[] newGrad = new double[dim];

      double alphaScaled = alpha / n; // L2 regularization

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

      // early exit vector
      double[] prevWeights = new double[dim];
      for (int j = 0; j "lt" dim; ++j)
        prevWeights[j] = this.weights[j];
      
      for (int epoch = 0; epoch "lt" maxEpochs; ++epoch)
      {
        Shuffle(indices, this.rnd);

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

          double err = predY - actualY;
          for (int j = 0; j "lt" dim; ++j) // each weight
          {
            newGrad[j] = x[j] * err;
            totGrad[j] = 
              totGrad[j] - gradMem[ii][j] + newGrad[j];
            gradMem[ii][j] = newGrad[j];

            this.weights[j] -= 
              lrnRate *  (totGrad[j] / n);
            this.weights[j] -= 
              lrnRate * alphaScaled * this.weights[j];
          }
          this.bias -= lrnRate * err;

        }
        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 early-stop after each epoch
        // if max_change_in_wts / max_weights "lt" 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 = weightDeltas[0];
        for (int j = 0; j "lt" dim; ++j)
          if (weightDeltas[j] "gt" maxChange)
            maxChange = weightDeltas[j];
        double maxWeight = this.weights[0];
        for (int j = 0; j "lt" dim; ++j)
          if (this.weights[j] "gt" maxWeight)
            maxWeight = this.weights[j];
        if (maxWeight != 0.0 &&
          (maxChange / maxWeight) "lt" stopTol)
        {
          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;
    }

    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
Posted in Machine Learning | Leave a comment

Implementing a StandardScaler for Machine Learning Regression Using C#

I rarely use L2 regularization for linear regression (“ridge regression”) and I never use L1 regularization for linear regression (“lasso regression”). But, in the rare situations when I use L2 ridge regression, it’s necessary (in theory) to scale the training and test X predictor data so that all model weights are treated equally by L2 weight reduction.

One morning I realized it had been many months since I used a C# standard scaler, so I implemented one. A standard scaler converts each column of training data into a new column where the values have a mean of 0.0 (“centered”) and variance of 1.0 (“unit variance”).

This is done using the computation xij’ = (xij – u) / sd, where xij’ is the scaled value at row i col j, xij is the unscaled value at row i col j, u is the mean of the column values, sd is the (population) standard deviation of the column values.

The output of my demo:

Begin StandardScaler with C# demo

Loading synthetic train (50) data

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

First three train y:
  0.4840
  0.1568
  0.8054

Creating standard scaler for trainX
Done

Column means:
 -0.0217  0.0201 -0.0283  0.0532  0.1676
Column variances:
  0.3665  0.3349  0.3905  0.2815  0.3713

Transforming data
Done

First three scaled train X:
 -0.2383  0.7266 -1.5546 -0.8454 -1.4346
  0.1640 -0.3140  0.6380 -1.2145  0.9661
 -1.5254  0.5543 -0.2194  0.1209 -1.4554

Unscaling scaled data

First three original 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

End demo

I used the scikit-learn library StandardScaler as a design guide for my C# class. This design stores column variances rather than column standard deviations — I’m not exactly sure why, but it’s not important.

The key calling statements are:

// read trainX matrix from file
StandardScaler scaler = new StandardScaler();
scaler.Fit(trainX); 
double[][] scaledTrainX = scaler.Transform(trainX);

In a non-demo scenario, you’d also scale the test data using the scaler that was fitted to the test data, so that there’s no data leakage.

Good fun.



Galaxy Science Fiction was published from 1950 to 1980. It was the leading science fiction magazine of its time. The magazine had a standardized cover format with beautiful art by a few dozen different artists. Here are three nice covers by artist Mel Hunter (1927-2004).


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

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

namespace Scaler
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin StandardScaler" +
        " with C# demo ");

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

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

      Console.WriteLine("\nCreating standard " +
        "scaler for trainX ");
      StandardScaler scaler = new StandardScaler();
      scaler.Fit(trainX);
      Console.WriteLine("Done ");

      Console.WriteLine("\nColumn means: ");
      VecShow(scaler.means, 4, 8);
      Console.WriteLine("Column variances: ");
      VecShow(scaler.variances, 4, 8);

      Console.WriteLine("\nTransforming data ");
      double[][] scaledTrainX = scaler.Transform(trainX);
      Console.WriteLine("Done ");

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

      Console.WriteLine("\nUnscaling scaled data ");
      double[][] origX = 
        scaler.InverseTransform(scaledTrainX);
     
      Console.WriteLine("\nFirst three original train X: ");
      for (int i = 0; i "lt" 3; ++i)
        VecShow(origX[i], 4, 8);

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

  } // Program

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

  public class StandardScaler
  {
    public double[] means;
    public double[] variances;

    public StandardScaler()
    {
      this.means = new double[0];
      this.variances = new double[0];
    }

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

    public void Fit(double[][] dataX)
    {
      int n = dataX.Length;
      int dim = dataX[0].Length;

      this.means = new double[dim];
      this.variances = new double[dim];

      for (int j = 0; j "lt" dim; ++j) // each col
      {
        double sum = 0.0;
        for (int i = 0; i "lt" n; ++i)
          sum += dataX[i][j];
        this.means[j] = sum / n;
      }

      for (int j = 0; j "lt" dim; ++j) // each col
      {
        double sum = 0.0;
        for (int i = 0; i "lt" n; ++i)
          sum += (dataX[i][j] - this.means[j]) *
            (dataX[i][j] - this.means[j]);
        this.variances[j] = sum / n;
      }
    }

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

    public double[][] Transform(double[][] dataX)
    {
      // x' = (x - u) / sd
      int n = dataX.Length;
      int dim = dataX[0].Length;

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

      for (int j = 0; j "lt" dim; ++j)
      {
        for (int i = 0; i "lt" n; ++i)
        {
          double x = dataX[i][j];
          double u = this.means[j];
          double sd = Math.Sqrt(this.variances[j]);
          result[i][j] = (x - u) / sd;
        }
      }
      return result;
    }

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

    public double[][] InverseTransform(double[][] scaledX)
    {
      // x = (x' * sd) + u

      int n = scaledX.Length;
      int dim = scaledX[0].Length;

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

      for (int j = 0; j "lt" dim; ++j)
      {
        for (int i = 0; i "lt" n; ++i)
        {
          double x = scaledX[i][j];
          double u = this.means[j];
          double sd = Math.Sqrt(this.variances[j]);
          result[i][j] = (x * sd) + u;
        }
      }
      return result;
    }

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

  } // class StandardScaler

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

} // ns

Test data:

# synthetic_train_50.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
Posted in Machine Learning | Leave a comment

Checking Machine Learning Training Data for Multicollinearity Using VIF (Variance Inflation Factor) from Scratch Python

In machine learning, if training data is multicollinear, the resulting model will likely be poor. 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 (well in weird scenarios, a VIF value could be less than one). Briefly, if all column VIF values are less than about 7.0, the data is probably OK.

if VIF is close to 1.0, the column is not correlated.
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 the 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.

I put together a demo using Python and the scikit library. 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 the demo program is:

Begin variance inflation factor demo

Loading synthetic (20) normal data

First two X:
[-0.1660  0.4406 -0.9998 -0.3953 -0.7065]
[ 0.0776 -0.1616  0.3704 -0.5911  0.7562]

Begin VIF analysis
col =   0 | vif = 1.1980
col =   1 | vif = 1.4591
col =   2 | vif = 1.2345
col =   3 | vif = 1.3025
col =   4 | vif = 1.2120

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

First two X:
[-0.1660  0.4406  0.1096 -0.3953 -0.7065]
[ 0.0776 -0.1616 -0.0045 -0.5911  0.7562]

Begin VIF analysis
col =   0 vif = 25546262.9389
col =   1 vif = 6023299.3886
col =   2 vif = 30951889.1370
col =   3 vif = 1.2937
col =   4 vif = 1.2117

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.

No moral to this blog post. Just 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 want a good relationship between the hero and the main actress.

I’m a huge fan of science fiction movies from the 1950s and 1960s. Here are posters of two films that were good, but they could have been great if the chemistry between the hero and the main lady were better.

Left: In “Crack in the World” (1965), scientists create a project to drill to the Earth’s magma center to gain a source of unlimited heat, and therefore unlimited energy. The plan involves firing a thermonuclear missile into a hole. This was not a good idea, to put it mildly. The chemistry between Dr. Rampion (actor Kieron Moore) and the wife of his boss, Dr. Sorensen (actress Janette Scott) was, well, one with no chemistry. But I give the movie a B grade anyway.

Right: In “The Day the Earth Caught Fire” (1961), The U.S. and the Soviets unknowingly explode nuclear test weapons at the same time on the same day. This was not a good idea, to put it mildly. The Earth is knocked out of orbit, towards the Sun. Only exploding every nuclear device on the planet simultaneously might save humanity. The chemistry between newspaper reporter Peter Stenning (actor Edward Judd) and office worker Jeannie Craig (actress Janet Munro — she’s one of my sci fi favorites) was awkward and unconvincing. But I give the movie a B- grade anyway.


Demo program:

# variance_inflation_factor.py

import numpy as np
from sklearn.linear_model import LinearRegression

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

def vif(data, i):
  # vif = 1.0 / (1.0 - R2) if col [i] is dependent variable
  X = np.delete(data, i, axis=1) # all cols except i
  y = data[:,i]

  model = LinearRegression()
  model.fit(X, y)
  r2 = model.score(X, y)
  result = 1.0 / (1.0 - r2)
  return result

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

print("\nBegin variance inflation factor demo ")

print("\nLoading synthetic (20) normal data ")
train_X = \
  np.loadtxt(".\\Data\\synthetic_train_20.txt",
  usecols=[0,1,2,3,4], delimiter=",")
 
print("\nFirst two X: ")
for i in range(2):
  print(train_X[i])

print("\nBegin VIF analysis ")

for c in range(len(train_X[0])):
  z = vif(train_X, c)
  print("col = %3d | vif = %0.4f " % (c, z))

print("\nLoading synthetic (20) mulicollinear data ")
print("(col[2] = 2.0 * col[0] + col[1] + rnd) ")
train_X = \
  np.loadtxt(".\\Data\\synthetic_train_20_collinear.txt",
  usecols=[0,1,2,3,4], delimiter=",")

print("\nFirst two X: ")
for i in range(2):
  print(train_X[i])

print("\nBegin VIF analysis ")

for c in range(len(train_X[0])):
  z = vif(train_X, c)
  print("col = %3d vif = %0.4f " % (c, z))
print("\nEnd demo ")

First, normal, 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, 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
Posted in Machine Learning, Scikit | Leave a comment