Gradient Boost Regression Using scikit Applied to the Diabetes Dataset – Poor Results As Expected

I write code almost every day. Like many skills, writing code is something that must be practiced, and anyway, I just enjoy writing code. One morning, I figured I’d run the well-known Diabetes Dataset through a from-scratch C# implementation of gradient boost regression. I realized that this was going to be a non-trivial exploration. So I decided to run the Diabetes Dataset through the scikit GradientBoostingRegressor module, so that I could get some baseline sanity-check results to see if my future scratch C# system is on the right track.

Based on previous experiments with linear regression, quadratic regression, neural network regression, kernel ridge regression, random forest regression, and AdaBoost regression, I was pretty sure the scikit GradientBoostingRegressor model would give poor prediction accuracy, and that’s what happened.

The raw Diabetes Dataset looks like:

59, 2, 32.1, 101.00, 157,  93.2, 38, 4.00, 4.8598, 87, 151
48, 1, 21.6,  87.00, 183, 103.2, 70, 3.00, 3.8918, 69,  75
72, 2, 30.5,  93.00, 156,  93.6, 41, 4.00, 4.6728, 85, 141
. . .

Each line represents a patient. The first 10 values on each line are predictors. The last value on each line is the target value (a diabetes metric) to predict. The predictors are: age, sex, body mass index, blood pressure, serum cholesterol, low-density lipoproteins, high-density lipoproteins, total cholesterol, triglycerides, blood sugar. There are 442 data items.

The sex encoding isn’t explained anywhere but I suspect male = 1, female = 2 because there are 235 1 values and 206 2 values).

Note that this Diabetes Dataset, which is included as an example dataset in the Python language scikit-learn library, is not the same as the Pima Diabetes Dataset from the UCI dataset repository.

For the sanity-check scikit version, I just use the built-in scikit dataset, with the model_selection train_test_split() function which will scale the predictors, and use 75% of the 442 items for training (331 items) and 25% for test (111 items). To keep the target y values in roughly the same range as the predictors, I divide the y values by 1000, but this doesn’t affect model accuracy.

The output of the demo is:

Begin scikit Gradient Boost demo

First three X predictors:
[[-0.0491 -0.0446 -0.0569 -0.0435 -0.0456 -0.0433  0.0008
  -0.0395 -0.0119  0.0155]

 [-0.0527 -0.0446 -0.0558 -0.0367  0.0892 -0.0032  0.0081
   0.0343  0.1324  0.0031]

 [-0.0927  0.0507 -0.0903 -0.0573 -0.0250 -0.0304 -0.0066
  -0.0026  0.0241  0.0031]]

First three y targets:
0.0680
0.1090
0.0940

Setting n_estimators = 100
Setting max_depth = 3
Setting learning rate = 0.1000

Training Gradient Boost model
Done

Accuracy train (within 0.10): 0.4079
Accuracy test (within 0.10): 0.1982

MSE train: 0.0008
MSE test: 0.0039

R2 train: 0.8785
R2 test: 0.2165

Predicting for x =
[-0.0491 -0.0446 -0.0569 -0.0435 -0.0456
 -0.0433  0.0008 -0.0395 -0.0119 0.0155]
Predicted y = 0.0961

End demo

These poor results were similar to the results I got using other regression models. I have done many experiments with the Diabetes Dataset and I’ve concluded the the default target value in the last column (a patient diabetes score) simply cannot be predicted well. But the variables in columns [4], [5], [6], [7], and [8] can be meaningfully predicted from the other columns.



Gradient boost regression was new and exciting when it was developed in 1999-2001. Twenty-five years later, it’s reaching “vintage machine learning technique” status.

I have always been a big fan of vintage space art from the 1950s. Here are a couple of examples by two of my favorite artists. Left: A lunar outpost by Alexander Leydenfrost (1888-1961). Right: An orbiting space station in operation by Alex Schomburg (1905-1998).


Demo program. Replace “lt” in the accuracy() function with the Boolean less-than operator. (My blog editor chokes on symbols).

# scikit_gradient_boost.py
# demo using Diabetes Dataset

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor

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

def accuracy(model, data_X, data_y, pct_close):
  # assumes model has a predict(X)
  n = len(data_X)
  n_correct = 0; n_wrong = 0
  for i in range(n):
    x = data_X[i].reshape(1,-1)  # make it a matrix
    y = data_y[i]
    y_pred = model.predict(x)[0]  # predict() expects 2D

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

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

def MSE(model, data_X, data_y):
  # scratch. could use sklearn.metrics.mean_squared_error()
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i].reshape(1,-1)
    y = data_y[i]
    y_pred = model.predict(x)[0]
    sum += (y - y_pred) * (y - y_pred)

  return sum / n

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

def main():
  print("\nBegin scikit Gradient Boost demo ")

  np.set_printoptions(precision=4, suppress=True,
    floatmode='fixed')
  np.random.seed(0)  # not used this version

  # use built-in scikit diabetes data
  X, y = load_diabetes(return_X_y=True, scaled=True)
  train_X, test_X, train_y, test_y = \
    train_test_split(X, y, random_state=0)  # 25% test

  # scale y to keep MSE values more interpretable
  train_y /= 1000
  test_y /= 1000

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

  n_ests = 100  # these are default values
  max_d = 3
  lrn_rate = 0.10

  print("\nSetting n_estimators = " + \
    str(n_ests))
  print("Setting max_depth = " + \
    str(max_d))
  print("Setting learning rate = %0.4f" % lrn_rate)

  # GradientBoostingRegressor(*, loss='squared_error',
  # learning_rate=0.1, n_estimators=100, subsample=1.0,
  # criterion='friedman_mse', min_samples_split=2,
  # min_samples_leaf=1, min_weight_fraction_leaf=0.0,
  # max_depth=3, min_impurity_decrease=0.0, init=None,
  # random_state=None, max_features=None, alpha=0.9,
  # verbose=0, max_leaf_nodes=None, warm_start=False,
  # validation_fraction=0.1, n_iter_no_change=None,
  # tol=0.0001, ccp_alpha=0.0)
  # SHEESH!

  model = \
    GradientBoostingRegressor(
    n_estimators=n_ests,
    max_depth=max_d,
    learning_rate=lrn_rate,
    random_state=0)

  print("\nTraining Gradient Boost model ")
  model.fit(train_X, train_y)
  print("Done " )

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

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

  r2_train = model.score(train_X, train_y)
  print("\nR2 train: %0.4f " % r2_train)
  r2_test = model.score(test_X, test_y)
  print("R2 test: %0.4f " % r2_test)

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

  print("\nEnd demo ")

if __name__ == "__main__":
  main()
Posted in Machine Learning, Scikit | Leave a comment

Linear Ridge Regression From Scratch Using Python With Closed Form Training

One morning before work, I noticed that it had been quite a few months since I last implemented linear ridge regression, from scratch, using Python, with closed form training instead of SGD training. I figured I’d do so, guided somewhat by the style of the scikit-learn Ridge module.

Ridge regression (sometimes called, and the same as, linear ridge regression, but is much different from kernel ridge regression) is just basic linear regression with L2 regularization. L2 regularization is applied during training and discourages model weights from becoming extremely large, because large model weights indicate likely model overfitting.

Mathematically, L2 regularization adds a penalty equal to the sum of the squared weights in the model. When you train a linear ridge regression model using SGD, the implementation is straightforward. But when you train using a closed form algorithm, typically left pseudo-inverse via normal equations with Cholesky inverse, or relaxed Moore-Penrose pseudo-inverse with SVD inverse, to implement L2 regularization, you add the L2 alpha constant to the diagonal elements of the pertinent matrix just before computing the pertinent inverse (Cholesky or SVD).

Adding the alpha L2 constant to the pertinent matrix does two things. First, it introduces L2 regularization (not obvious at all, but mathematically correct). Second, it conditions the matrix so the Cholesky inverse will not fail.


Some Minor Details

If you use SGD training, you modify the gradient term for each weight just before the update statement. Simple, but SGD requires a learning rate and a max epochs, which must be determined by trial and error.

In math literature, the L2 constant is called lambda, but in computer science we use alpha to avoid colliding with the lambda keyword that is found in many programming languages.

There is also L1 regularization, but I never use it — explaining why is a non-trial topic.

In math notation, the number of rows in a matrix is usually called m and the number of columns is n. I use non-standard notation of n for the number of rows, and dim for the number of columns.

My demo implementation uses explicit for-loops for clarity. Numpy supports all kinds of bulk operations on vectors and matrices, which is more efficient, but much less clear.

Years ago, the error term was sometimes computed as pred_y – target_y and sometimes computed as target_y – pred_y. The pred_y – target_y form seems universal now. If you use target_y – pred_y, then the -= in the weight update becomes += (which always seemed a bit more logical to me, but I sadly was never elected King of All Machine Learning).

The scikit Ridge module uses a clever mini-algorithm to automatically estimate a good value for the learning rate, but I suspect that the clever mini-algorithm algorithm works well only with the weird SAG training algorithm. You can always use a simple grid search to find a good learning rate.

I use A[i][j] matrix syntax rather than the slightly more efficient A[i,j] syntax — just a personal preference.

To evaluate the from-scratch ridge regression model and the scikit Ridge model, I implemented a program-defined mse() function and an accuracy() function. I allowed for the fact that all scikit model predict() methods require a matrix as input and return a vector as output, rather than accepting a vector as input and returning a single value as output.


Mathematically, the bias and the weights are computed directly by bias_wts = inv(Xt * X) * Xt * train_y where X is the design matrix derived from the train_X matrix, Xt is the transpose of X, and * is matrix multiplication or matrix-to-vector multiplication.

The key training code looks like:

  # compute design matrix from train X
  col = np.ones((n,1))
  design_X = np.hstack((col, train_X))

  # compute Xt * X
  Xt = np.transpose(design_X)
  XtX = np.matmul(Xt, design_X)

  # condition Xt*X and add L2 before calling inv()
  for i in range(dim):
    XtX[i][i] += alpha

  inv_XtX = self.cholesky_inv(XtX)
  tmp = np.matmul(inv_XtX, Xt)
  bias_weights = np.dot(tmp, train_y)

  # extract to model
  self.bias = bias_weights[0];
  for i in range(1,len(bias_weights)):
    self.weights[i-1] = bias_weights[i]

Somewhat strangely, there is no-built-in Cholesky inverse function, so I implemented one using the built-in Cholesky decomposition function. You can use any general inverse function, but a conditioned Xt * X matrix is square, symmetric, positive definite, which is a special case which can be easily inverted using Cholesky.

For my demo, I used a set of synthetic data that 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 first five values on each line are predictors. The last value on each line is the target y value to predict. The data has been pre-scaled so that all predictors are between -1 and +1. When using ridge regression, you should scale data so that all weights are penalized equally. The scikit StandardScaler module can do this if your data is not pre-scaled.

There are 200 training items and 40 test items. The data has complex non-linear relations so linear regression, ridge or plain, will not work very well.

The key parts of the demo output are:

Begin ridge regression scratch Python
Loading train (200) and test (40) data

Creating scratch Python ridge regression model
Done

Training model with left pinv via normal eqs
Setting L2 alpha = 1.000000
Done.

Model weights:
[-0.2617  0.0330 -0.0452  0.0350 -0.1145]
Model bias = 0.3600

Evaluating model

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

MSE train = 0.00259
MSE test = 0.00192

To sort-of validate my demo system, I ran the data through the scikit Ridge module. The model produced by scikit was almost identical to the from-scratch model:

Model weights:
[-0.2618  0.0331 -0.0453  0.0353 -0.1132]
Model bias = 0.3618

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

MSE train = 0.00259
MSE test = 0.00194

One weird thing about the scikit Ridge module is that it doesn’t use alpha directly. Instead it scales alpha by dividing it by the number of training items. This is not documented and it’s not clear to me exactly what scikit Ridge is doing.

The scikit Ridge module can use one of seven training algorithms: ‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga, ‘lbfgs’. The ‘auto’ means use some mini-algorithm to guess the best training algorithm. My demo specified ‘cholesky’ so the model is as close to the from-scratch version as possible.



Machine learning linear regresssion techniques, including ridge regression, were developed in the late 1950s and early 1960s. The space race to the moon was also going on in those years. It was a thrilling time to be alive. And I love science fiction movies from that era.

Two of my favorite science fiction movies from the 1950s feature alien spacecraft that have a spherical design.

Left: In “The Man from Planet X” (1951), an alien spaceship lands in the remote Scottish moors. The alien on board is a scouting party for an invasion from an approaching planet X. The alien uses mind control to keep local villagers at bay, but ultimately the alien is destroyed by the UK military and the invasion is thwarted. Objectively this movie isn’t very good, but I give it my personal B grade. Some of the scenes on the moors absolutely terrified me when I was young.

Right: In “It Came from Outer Space” (1957), an astonomer and his girlfriend observe a huge spherical object crash into the desert and bury itself underground. Soon local townspeople start disapperaing, then reappearing and acting strangely. It turns out that the aliens are benign and just want to repair their ship. The story ends well for everyone. Considered a classic. I give this movie my personal A grade.


Demo program. Replace the “lt” in the accuracy() function with the Boolean less-than operator.

# ridge_regression_closed_form_scratch.py
# left pseudo-inverse via normal equations training

import numpy as np

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

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

    if scikit == False:
      pred_y = model.predict(x)
    else:
      pred_y = model.predict([x])[0]

    if np.abs(pred_y - 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, scikit=False):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i]
    y = data_y[i]
    if scikit == False:
      pred_y = model.predict(x)
    else:
      pred_y = model.predict([x])[0]
    diff = pred_y - y
    sum += diff * diff
  return sum /n

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

class RidgeRegressor:
  # ridge regression using SGD
  def __init__(self, seed=0):
    self.weights = None
    self.bias = 0.0
    self.rnd = np.random.RandomState(seed) # not used

  def predict(self, x):
    # x is a vector
    sum = 0.0
    for j in range(len(x)):
      sum += self.weights[j] * x[j]
    sum += self.bias
    return sum

  def train(self, train_X, train_y, alpha=0.01):
    # bias_wts = inv(Xt * X) * Xt * y
    n = len(train_X)
    dim = len(train_X[0])  # aka m,n in math
    self.weights = np.zeros(dim)

    # compute design matrix from X
    col = np.ones((n,1))
    design_X = np.hstack((col, train_X))

    # compute Xt * X
    Xt = np.transpose(design_X)
    XtX = np.matmul(Xt, design_X)

    # condition Xt*X before calling inv()
    for i in range(dim):
      XtX[i][i] += alpha

    # inv_XtX = np.linalg.inv(XtX) # no direct Cholesky
    inv_XtX = self.cholesky_inv(XtX)
    tmp = np.matmul(inv_XtX, Xt)
    bias_weights = np.dot(tmp, train_y)

    # extract to model
    self.bias = bias_weights[0];
    for i in range(1,len(bias_weights)):
      self.weights[i-1] = bias_weights[i]
    return  # all done

  def cholesky_inv(self, A):
    # weirdly, there is no direct Cholesky inverse
    L = np.linalg.cholesky(A) # Cholesky decomp
    Linv = np.linalg.inv(L) # sadly, not optimized
    Ainv = np.matmul(Linv.T, Linv)
    return Ainv

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

print("\nBegin ridge regression scratch Python ")

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

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

print("\n========================== ")

print("\nCreating scratch Python ridge regression model ")
model = RidgeRegressor()
print("Done ")

print("\nTraining model with left pinv via normal eqs ")
alpha = 1.0
print("Setting L2 alpha = %0.6f " % alpha)
model.train(train_X, train_y, alpha)
print("Done. " )

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

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.5f " % mse_train)
print("MSE test = %0.5f " % mse_test)

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

print("\n========================== ")

from sklearn.linear_model import Ridge

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

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

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

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

acc_train = \
  accuracy(model, train_X, train_y, 0.10, scikit=True)
acc_test = \
  accuracy(model, test_X, test_y, 0.10, scikit=True)
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, scikit=True)
mse_test = mse(model, test_X, test_y, scikit=True)
print("\nMSE train = %0.5f " % mse_train)
print("MSE test = %0.5f " % mse_test)

print("\nEnd ridge regression 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

Right Pseudo-Inverse via Normal Equations Implemented Using C#

Bottom line: Right pseudo-inverse is virtually never used in machine learning scenarios, but I put together a demo anyway.

One morning before work, I was waiting for the rain to stop so I could give my two dogs their early walk. I figured I’d put together a demo of computing the right pseudo-inverse of a matrix using the C# language. I have implemented left-pseudo inverse many times — it’s a useful technique for training linear models — but I haven’t used right-pseudo-inverse in recent memory.

Briefly, left-pseudo-inverse is helpful in machine learning scenarios where a matrix of training data has more rows than columns. The right-pseudo inverse is helpful in scientific scenarios where a matrix of some sort of data has more columns than rows.

If you have a square matrix A, the inverse of A, Ainv, is a matrix so that A * Ainv = Ainv * A = I, where * is matrix multiplication and I is the identity matrix (1.0s on the upper-left to lower-right diagonal, 0.0s elsewhere). Not all square matrices have an inverse.

For non-square matrices, it’s useful to compute a pseudo-inverse. In machine learning scenarios, the matrix to invert is usually training data, or a design matrix (training data with a leading column of 1.0s added). In these scenarios, the number of rows (number of training items) is almost always greater than the number of columns (predictor variables).

There are two main categories of pseudo-inverse algorithms used in machine learning training scenarios: 1.) relaxed Moore-Penrose pseudo-inverse (using either QR or SVD inverse, each with several variations), 2.) left pseudo-inverse (almost always using Cholesky inverse). Relaxed Moore-Penrose pseudo-inverse is extremely complex. Left pseudo-inverse is simpler but works only when nRows is greater than nCols — which is the case for machine learning scenarios.

If you have a matrix A with more columns than rows (sometimes occurs in scientific programming but rarely in machine learning scenarios), the right pseudo-inverse via normal equations is computed as right-pinv(A) = At * inv(A * At). Here inv() is any regular matrix inverse function, but because A * At is square symmetric positive definite (assuming you condition the diagonal by adding a small constant, typically about .0e-8), you can apply the relatively simple Cholesky inverse function.

In code:

public class Cholesky
{
  public static double[][] MatRightPseudoInv(double[][] A)
  {
    // right-pinv(A) = At * inv(A * At)
    double[][] At = MatTranspose(A);
    double[][] AAt = MatProduct(A, At);
    for (int i = 0; i < A.Length; ++i) // condition AAt
      A[i][i] += 1.0e-8;
    double[][] AAtinv = MatInvCholesky(AAt);
    double[][] result = MatProduct(At, AAtinv);
    return result;
  }
  . . . 

I wrote a demo program that generates random matrices, computes the right pseudo-inverse, and checks the result. The right pseudo-inverse worked reasonably well (precision to within 1.0e-6 but no better than that).

The main problem with both left-pseudo inverse and right pseudo-inverse is the computation of At * A or A * At. If there are many small values, or many large values in At, then A * At can generate arithmetic overflow or underflow.



The 1930s, 40s, and 50s were the golden age of science fiction magazines. Big bugs made frequent cover appearances.

Left: Amazing Stories, May 1941.
Center: Thrilling Wonder Stories, December 1939.
Right: Fantastic Adventures, February 1940.


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

using System.IO;

namespace MatrixRightPseudoInverse
{
  internal class MatrixRightPseudoInverseProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin right pseudo-inverse" +
        " via normal equations (Cholesky) ");
      Console.WriteLine("Used when nCols gt nRows (rare) ");

      int minRows = 2; int maxRows = 10;
      int minCols = 10; int maxCols = 1000;
      Random rnd = new Random(0);
      int nTrials = 1000;
      int nPass = 0; int nFail = 0;
      for (int i = 0; i "lt" nTrials; ++i)
      {
        Console.WriteLine("\n==========");
        Console.WriteLine("trial " + i);
        double[][] A =
          MakeRandomMatrix(minRows, maxRows,
          minCols, maxCols, rnd);

        if (A.Length "gt" A[0].Length)
          continue;  // right pseudo-inv doesn't apply
        double[][] Apinv = Cholesky.MatRightPseudoInv(A);

        // check (A * Apinv) * A = A
        double[][] AApinv = MatProduct(A, Apinv);
        double[][] C = MatProduct(AApinv, A);

        //Console.WriteLine("\nA = ");
        //MatShow(A, 4, 9);
        //Console.WriteLine("\nApinv = ");
        //MatShow(Apinv, 4, 9);
        //Console.WriteLine("\n(A * Apinv) * A = ");
        //MatShow(C, 4, 9);

        if (MatAreEqual(C, A, 1.0e-6) == true)
        {
          Console.WriteLine("pass");
          ++nPass;
        }
        else
        {
          Console.WriteLine("FAIL");
          Console.WriteLine("nRows = " + A.Length +
            " nCols = " + A[0].Length);
          Console.ReadLine();
          ++nFail;
        }
        Console.WriteLine("==========");
     } // nTrials

      Console.WriteLine("\nNumber pass = " + nPass);
      Console.WriteLine("Number fail = " + nFail);

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

    // ------------------------------------------------------
    // helpers: MakeRandomMatrix(), MatProduct(), 
    // MatShow(), MatAreEqual()
    // ------------------------------------------------------

    public static double[][] MakeRandomMatrix(int minRows,
      int maxRows, int minCols, int maxCols, Random rnd)
    {
      int nRows = rnd.Next(minRows, maxRows);  // [2,maxN)
      int nCols = rnd.Next(minCols, maxCols);
      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols];

      double lo = -10.0; double hi = 10.0;
      for (int r = 0; r "lt" nRows; ++r)
        for (int c = 0; c "lt" nCols; ++c)
          result[r][c] = (hi - lo) * rnd.NextDouble() + lo;
      return result;
    }

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

    public static double[][] MatProduct(double[][] A,
      double[][] B)
    {
      int aRows = A.Length; int aCols = A[0].Length;
      int bRows = B.Length; int bCols = B[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

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

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

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

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

    public static bool MatAreEqual(double[][] A, double[][] B,
      double epsilon)
    {
      int nRows = A.Length; int nCols = A[0].Length;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          if (Math.Abs(A[i][j] - B[i][j]) "gt" epsilon)
            return false;
      return true;
    }

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

  } // Program

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

  public class Cholesky
  {
    public static double[][] MatRightPseudoInv(double[][] A)
    {
      // right-pinv(A) = At * inv(A * At)
      double[][] At = MatTranspose(A);
      double[][] AAt = MatProduct(A, At);
      // condition AAt
      for (int i = 0; i "lt" A.Length; ++i)
        A[i][i] += 1.0e-8;
      double[][] AAtinv = MatInvCholesky(AAt);
      double[][] result = MatProduct(At, AAtinv);
      return result;
    }
    public static double[][] MatLeftPseudoInv(double[][] A)
    {
      // left-pinv(A) = inv(At * A) * At
      double[][] At = MatTranspose(A);
      double[][] AtA = MatProduct(At, A);
      // condition AtA
      for (int i = 0; i "lt" A.Length; ++i)
        A[i][i] += 1.0e-8;
      double[][] AtAinv = MatInvCholesky(AtA);
      double[][] result = MatProduct(AtAinv, At);
      return result;
    }

    private static double[][] MatInvCholesky(double[][] A)
    {
      // A must be square symmetric positive definite
      // A = L * Lt, inv(A) = inv(Lt) * inv(L)
      double[][] L = MatDecompCholesky(A); // lower tri
      double[][] Lt = MatTranspose(L); // upper tri
      double[][] invL = MatInvLowerTri(L);
      double[][] invLt = MatInvUpperTri(Lt);
      double[][] result = MatProduct(invLt, invL);
      return result;
    }

    private static double[][] MatDecompCholesky(double[][] A)
    {
      // A must be square, symmetric, positive definite
      int m = A.Length; int n = A[0].Length;  // m == n
      double[][] L = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        L[i] = new double[n];

      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lte" i; ++j)
        {
          double sum = 0.0;
          for (int k = 0; k "lt" j; ++k)
            sum += L[i][k] * L[j][k];
          if (i == j)
          {
            double tmp = A[i][i] - sum;
            if (tmp "lt" 0.0)
              throw new
                Exception("decomp Cholesky fatal");
            L[i][j] = Math.Sqrt(tmp);
          }
          else
          {
            if (L[j][j] == 0.0)
              throw new
                Exception("decomp Cholesky fatal ");
            L[i][j] = (A[i][j] - sum) / L[j][j];
          }
        } // j
      } // i

      return L;
    }

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

    private static double[][] MatInvLowerTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[k][j] -= result[i][j] * A[k][i];
          }
          result[k][j] /= (A[k][k] + 1.0e-8);
        }
      }
      return result;
    }

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

    private static double[][] MatInvUpperTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[j][k] -= result[j][i] * A[i][k];
          }
          result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
        }
      }
      return result;
    }

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

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

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

    private static double[][] MatIdentity(int n)
    {
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      return result;
    }

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

    private static double[][] MatTranspose(double[][] M)
    {
      int nRows = M.Length; int nCols = M[0].Length;
      double[][] result = MatMake(nCols, nRows);  // note
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[j][i] = M[i][j];
      return result;
    }

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

    private static double[][] MatProduct(double[][] A,
      double[][] B)
    {
      int aRows = A.Length; int aCols = A[0].Length;
      int bRows = B.Length; int bCols = B[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

      double[][] result = MatMake(aRows, bCols);

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

  } // class Cholesky

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

} // ns
Posted in Machine Learning | Leave a comment

My Favorite Disney Animated Film From Each Decade

Disney was a huge part of my life. Growing up in Anaheim, in the late 1950s and early 1960s (shortly after Disneyland opened), the Matterhorn was visible from my house, and the Park was an immediate influence. In those days, admission to Disneyland was only $1.60 for adults and 60 cents for children, so we went frequently. Then, several years later, as a college student at UC Irvine in the 1970s, I worked at Disneyland in the evenings and weekends.

Here are my picks for my favorite Disney animated film in each decade, starting with the 1930s. I don’t include Pixar films.


1. Snow White and the Seven Dwarfs (1937) – This is an easy pick because it was the only animated film released in the 1930s. That said, it’s a brilliant film. It was the world’s first full-length animated movie and a bet-the-company gamble. If you can’t list the names of all seven dwarfs, shame on you.


2. Pinocchio (1940) – This film had absolutely ground-breaking technology. And it was quite scary too, with Monstro the evil whale and Pleasure Island where boys were turned into donkeys. A brilliant achievement. Honorable mentions: “Dumbo” (1941), “Bambi” (1942).


3. Sleeping Beauty (1959) – This film has my favorite art design and the single best female villain of all time: Maleficent. The dragon scene with green fire is a landmark in animation. Honorable mentions: “Cinderella” (1950), “Alice in Wonderland” (1951), “Peter Pan” (1953), “Lady and the Tramp” (1955) — all fantastic. The 1950s is the best decade for Disney animation in the Walt era prior to his death in 1966.


4. One Hundred and One Dalmatians (1961) – A wonderful story, even if the idea of raising Dalmatian puppies for their hides is rather grim. The highlight of this movie is the villainess, Cruella de Vil. Not a strong decade for Disney animation because Walt was focused on the “Mary Poppins” live-action movie and other projects. No honorable mentions for the 1960s.


5. The Rescuers (1977) – It took a while for Disney to regain momentum after the death of Walt in 1966. This movie isn’t highly regarded by many Disney experts but I give it a solid A grade. Mice Bernard and Bianca of the Rescue Aid Society save a young girl from Madame Medusa in the Louisiana swaps. “The Rescuers Down Under” (1990) was also excellent. No honorable mentions for the 1970s.


6. The Great Mouse Detective (1986) – Basil the detective (a Sherlock Holmesian mouse) and his friend Dr. Dawson, help a young mouse Olivia save Olivia’s father from the evil Professor Ratigan in Victorian England. The success of this film saved Disney animation, which was on the verge of bankruptcy after “The Dark Cauldron” (1985), which was a huge financial failure. Honorable mention: “The Little Mermaid” (1989), which most people other than me rate as far better than “Mouse Detective”.


7. Hercules (1997) – This film has a darker feel and more adult themes than many other Disney animated movies. Although this movie was successful, it wasn’t as successful as previous Disney blockbuster animated films of the 1990s: “Beauty and the Beast” (1991), “Aladdin” (1992), “The Lion King” (1994), and to a lesser degree “The Hunchback of Notre Dame” (1996). These four films all get honorable mention for the 1990s. Also highly rated by many people are “Mulan” (1998) and “Tarzan” (1999). An incredible decade for Disney animation.


8. Atlantis the Lost Empire (2001) – The story takes place in 1914 and tells the story of a scientific expedition that discovers Atlantis in spite of the evil mercenary Rourke. The movie only broke even at the box office — a major disappointment for Disney. There are several theories about why the movie wasn’t a success. All of the good animated films of this decade were done by Pixar — “The Incredibles”, “Monsters, Inc.”, “Finding Nemo”, “Cars”, “WALL-E”. A sad decade for Disney animation.


9. Zootopia (2016) – I avoided watching this movie for a long time because the basic idea sounded kind of lame — in a world of animals, a rookie police officer rabbit (Judy Hopps) and a con artist fox (Nick Wilde) work together to uncover a conspiracy involving the disappearance of predator animals. But when I finally watched, I was really impressed. I’d say this movie exceeds the sum of its parts. Unlike most people, I did not like “Frozen” (2013), or “Moana” (2016).


For the decade of the 2020s, so far (until mid-2026 as I write this post), there are no Disney animated films I rate as even better than so-so, led by the gosh-awful “Strange World” (2022) which lost $200 million, making it one of the biggest box-office bombs of all time. Others:

“Raya and the Last Dragon” (2021) – Not bad, but weak story and forgettable characters. And I didn’t like the animation style. This film lost money for the Studio. Grade = C.

“Encanto” (2021) – Again, not bad, but the story and characters are completely forgettable. Can you name any? This film lost money for the Studio. Grade = C.

“Wish” (2023) – Again, no real story to speak of. “A movie assembled by a focus group”. This film lost over $100 million for the Studio. Grade = C-.

“Moana 2” (2024) – Originally intended as a TV series. Animation outsourced, causing nearly 1,000 lost animator jobs. No good villain, no good songs, no good story. Grade = C.

“Mufasa: The Lion King” (2024) – Boring plot — not everything needs a backstory. Uninspired music. I really disliked the photorealistic animation. Creepy. Grade = C.

“Elio” (2025) – Weird. A giant box office bomb, and rightfully so. Lost about $40 million.

Disney animation has lost its way and has forgotten about producing wonderful stories that entertain, rather than the Studio’s current fixation on producing movies intended to influence social norms and generate merchandising opportunities. Shame on you Disney Studios. Right now you’re a disgrace to the memory of Walt Disney. Please regain your greatness and drop your self-destructive wokeness that alienates your audiences and loses huge amounts of money for your stockholders.


Posted in Top Ten | Leave a comment

“Quadratic Regression with Pseudo-Inverse Training Using C#” in Visual Studio Magazine

I wrote an article titled “Quadratic Regression with Pseudo-Inverse Training Using C#” in the May 2026 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2026/05/01/quadratic-regression-with-pseudo-inverse-training-using-csharp.aspx.

The goal of a machine learning regression problem is to predict a single numeric value. For example, you might want to predict a person’s credit rating based on age, annual salary, bank account balance, and so on. There are approximately a dozen common regression techniques. The most basic technique is called linear regression, or sometimes multiple linear regression, where the “multiple” indicates two or more predictor variables.

The form of a basic linear regression prediction model is y’ = (w0 * x0) + (w1 * x1) + . . + (wn * xn) + b, where y’ is the predicted value, the xi are n predictor values, the wi are model weights, and b is the bias.

Quadratic regression extends linear regression. The form of a quadratic regression model is y’ = (w0 * x0) + . . + (wn * xn) + (wj * x0 * x0) + . . + (wk * x0 * x1) + . . . + b. There are derived predictors that are the square of each original predictor, and interaction terms that are the multiplication product of all possible pairs of original predictors.

Compared to basic linear regression, quadratic regression can handle more complex data. Compared to the most powerful regression techniques such as gradient boosting regression, quadratic regression often has slightly worse prediction accuracy, but has much better model interpretability.

The output of the demo program presented in the article is:

Begin C# quadratic regression with pseudo-inverse
 (QR-Householder) training

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

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

First three train y:
  0.4840
  0.1568
  0.8054

Creating quadratic regression model

Starting pseudo-inverse training
Done

Model base weights:
 -0.2630  0.0354 -0.0421  0.0341 -0.1124

Model quadratic weights:
  0.0655  0.0194  0.0051  0.0047  0.0243

Model interaction weights:
  0.0043  0.0249  0.0071  0.1081 -0.0012 -0.0093
  0.0362  0.0085 -0.0568  0.0016

Model bias/intercept:   0.3220

Evaluating model
Accuracy train (within 0.10) = 0.8850
Accuracy test (within 0.10) = 0.9250

MSE train = 0.0003
MSE test = 0.0005

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

Predicted y = 0.4843

End demo

The model accuracy is quite good compared to many other regression techniques applied to the synthetic dataset. The model achieves 88.50% accuracy on the training data (177 out of 200 correct) and 92.50% accuracy on the test data (37 out of 40 correct). A prediction is scored as correct if it’s within 10% of the true target value.

Behind the scenes, the demo program trains the quadratic regression model using a technique called relaxed Moore-Penrose pseudo-inverse with QR inverse via the Householder algorithm. This is just one of many possible training algorithms. The different algorithms have different pros and cons and no algorithm works the best all the time.

Quadratic regression has a nice balance of prediction power and interpretability. The model weights/coefficients are easy to interpret. If the predictor values have been normalized to the same scale, larger magnitudes mean larger effect, and the sign of the weights indicate the direction of the effect.

Quadratic regression is not always effective — if it were, it would be used far more often than it is. Compared to basic linear regression, quadratic regression can sometimes provide a big improvement in model prediction accuracy for a relatively small investment in effort, and so it’s usually worth exploring.



Quadratic regression was developed in the 1950s and 1960s, about the same time as early space exploration. Galaxy Science Fiction was published from 1950 to 1980. It was the leading science fiction magazine of its time. I really enjoy looking at older covers to get a sense of the optimism and excitement in those days. Here are three nice 1952 covers by artist Jack Coggins (1911-2006).


Posted in Machine Learning | Leave a comment

AdaBoost Regression Using C# Applied to the Diabetes Dataset – Poor Results As Expected

I write code almost every day. Like many skills, writing code is something that must be practiced, and anyway, I just enjoy writing code. One morning before work, I figured I’d run the well-known Diabetes Dataset through an AdaBoost regression model. Based on previous experiments with linear regression, quadratic regression, neural network regression, kernel ridge regression, and random forest regression, I was pretty sure the C# Adaboost regression model would give poor prediction accuracy. And that’s what happened.

The raw Diabetes Dataset looks like:

59, 2, 32.1, 101.00, 157,  93.2, 38, 4.00, 4.8598, 87, 151
48, 1, 21.6,  87.00, 183, 103.2, 70, 3.00, 3.8918, 69,  75
72, 2, 30.5,  93.00, 156,  93.6, 41, 4.00, 4.6728, 85, 141
. . .

Each line represents a patient. The first 10 values on each line are predictors. The last value on each line is the target value (a diabetes metric) to predict. The predictors are: age, sex, body mass index, blood pressure, serum cholesterol, low-density lipoproteins, high-density lipoproteins, total cholesterol, triglycerides, blood sugar. There are 442 data items.

The sex encoding isn’t explained anywhere but I suspect male = 1, female = 2 because there are 235 1 values and 206 2 values).

Note that this Diabetes Dataset, which is included as an example dataset in the Python language scikit-learn library, is not the same as the Pima Diabetes Dataset from the UCI dataset repository.

I converted the sex values from 1,2 into 0,1. Then I applied divide-by-constant normalization by dividing the 10 predictor columns by (100, 1, 100, 1000, 1000, 1000, 100, 10, 10, 1000) and the target y values by 1000. The resulting encoded and normalized data looks like:

0.5900, 1.0000, 0.3210, . . . 0.1510
0.4800, 0.0000, 0.2160, . . . 0.0750
0.7200, 1.0000, 0.3050, . . . 0.1410
. . .

Normalization isn’t necessary for ensemble models like AdaBoost that use decision trees, but normalization doesn’t hurt. I split the 442-items into a 342-item training set and a 100-item test set.

AdaBoost uses a collection of decision trees. Each tree is trained using data that was poorly predicted by previous trees, so in theory, each tree gets better. The final predicted value is a weighted median of all the tree predictions.

I implemented an AdaBoost regression model using C#. The algorithm is known as AdaBoost.R2 (“adaptive boosting for regression, version 2”). I used 100 decision trees. I set the maximum depth of each tree to 5, the minimum samples parameter to 2 (standard default), and the minimum leaf parameter to 1 (standard default).

The output of the demo is:


Begin AdaBoost.R2 regression on Diabetes Dataset

Loading diabetes train (342) and test (100) normalized data
Done

First three train X:
  0.5900  1.0000  0.3210  . . .  0.4860  0.0870
  0.4800  0.0000  0.2160  . . .  0.3892  0.0690
  0.7200  1.0000  0.3050  . . .  0.4673  0.0850

First three train y:
  0.1510
  0.0750
  0.1410

Setting maxLearners = 100
Setting tree maxDepth = 5
Setting tree minSamples = 2
Setting tree minLeaf = 1

Training AdaBoost.R2 model
Done
Created 100 learners

Accuracy train (within 0.10): 0.2807
Accuracy test (within 0.10): 0.2065

MSE train: 0.0015
MSE test: 0.0035

Predicting for x =
  0.5900  1.0000  0.3210  . . .  0.4860  0.0870
Predicted y = 0.1779

End demo

These poor results were essentially the same as the results I got using the scikit AdaBoostRegressor module (there was slight differences because AdaBoost has a random component when selecting training items for each tree).

I have done many experiments with the Diabetes Dataset and I’ve concluded the the default target value in the last column (a patient diabetes score) simply cannot be predicted well. But the variables in columns [4], [5], [6], [7], and [8] can be meaningfully predicted from the other columns.



I have a fascination with machine learning. And I have a strange fascination with “chittering”. Several years ago, I noticed that when I watch movies with closed captioning enabled, there are an amazing number of scenes where an animal is making noises and the closed captioning reads “chittering”. I have seen this dozens of times. Here’s an example from “Superhero Movie” (2008), a very funny (to me) spoof of all the annoying superhero movies. On the left, the Dragonfly hero sees a crowd of people gathered. He rushes over and sees the crowd being fascinated by a monkey — while a super villain is robbing a bank and is being completely ignored.


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

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

namespace AdaBoostRegression
{
  internal class AdaBoostRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin AdaBoost.R2" +
        " regression on Diabetes Dataset ");

      // 1. load data
      Console.WriteLine("\nLoading diabetes train" +
        " (342) and test (100) normalized data");
      string trainFile = "..\\..\\..\\Data\\" +
        "diabetes_norm_train_342.txt";
      int[] colsX = new int[] { 0, 1, 2, 3, 4,
        5, 6, 7, 8, 9};
      int colY = 10;
      double[][] trainX =
        MatLoad(trainFile, colsX, ',', "#");
      double[] trainY =
        MatToVec(MatLoad(trainFile,
        new int[] { colY }, ',', "#"));

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

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

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

      // 2. create and train model
      int maxLearners = 100;
      int maxDepth = 5;
      int minSamples = 2;
      int minLeaf = 1;

      Console.WriteLine("\nSetting maxLearners = " +
        maxLearners);
      Console.WriteLine("Setting tree maxDepth = " +
        maxDepth);
      Console.WriteLine("Setting tree minSamples = " +
        minSamples);
      Console.WriteLine("Setting tree minLeaf = " +
        minLeaf);

      Console.WriteLine("\nTraining AdaBoost.R2 model ");
      AdaBoostRegressor model =
        new AdaBoostRegressor(maxLearners, maxDepth,
        minSamples, minLeaf, seed: 0); // master seed
      model.Train(trainX, trainY);
      Console.WriteLine("Done ");
      Console.WriteLine("Created " +
        model.learners.Count() + " learners ");

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

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

      // 4. use model to make a prediction
      double[] x = trainX[0];
      Console.WriteLine("\nPredicting for x = ");
      VecShow(x, 4, 8);
      double yPred = model.Predict(x);
      Console.WriteLine("Predicted y = " +
        yPred.ToString("F4"));
      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

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

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

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

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

  } // class Program

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

  class AdaBoostRegressor
  {
    public int maxLearners;
    public int maxDepth;
    public int minSamples;
    public int minLeaf;
    private Random rnd;
    private int seed; // master also passed to each tree
    public List"lt"DecisionTreeRegressor"gt" learners;
    public List"lt"double"gt" betas;

    public AdaBoostRegressor(int maxLearners, int maxDepth,
      int minSamples, int minLeaf, int seed)
    {
      this.maxLearners = maxLearners;
      this.maxDepth = maxDepth;
      this.minSamples = minSamples;
      this.minLeaf = minLeaf;
      this.rnd = new Random(seed);
      this.seed = seed;
      this.learners = new List"lt"DecisionTreeRegressor"gt"();
      this.betas = new List"lt"double"gt"();
      // each learner has a "model weight" derived
      // from its beta-confidence weight
    } // ctor

    public void Train(double[][] trainX, double[] trainY)
    {
      int N = trainX.Length;
      double[] weights = new double[N];
      for (int i = 0; i "lt" N; ++i)
        weights[i] = 1.0 / N;

      for (int t = 0; t "lt" this.maxLearners; ++t)
      {
        // 1. use weights to select a sample
        int[] sampleIdxs = this.SelectIdxs(weights, N);

        double[][] sampleX = ExtractRows(trainX,
          sampleIdxs);
        double[] sampleY = ExtractVals(trainY,
          sampleIdxs);

        // 2. construct regression machine t
        DecisionTreeRegressor lrnr =
          new DecisionTreeRegressor(this.maxDepth,
          this.minSamples, this.minLeaf, seed:this.seed);
        lrnr.Train(sampleX, sampleY);

        // 3. compute predicteds for all training items
        double[] predY = new double[N];
        for (int i = 0; i "lt" N; ++i)
          predY[i] = lrnr.Predict(trainX[i]);

        // 4. calculate loss each sample
        double[] absDiffs = new double[N];
        for (int i = 0; i "lt" N; ++i)
        {
          absDiffs[i] = Math.Abs(predY[i] - trainY[i]);
        }

        double D = 0.0;  // largest abs diff
        for (int i = 0; i "lt" N; ++i)
        {
          if (absDiffs[i] "gt" D)
            D = absDiffs[i];
        }

        double[] L = new double[N];
        for (int i = 0; i "lt" N; ++i)
          L[i] = absDiffs[i] / D;  // normalized

        // 5. calculate average loss
        double Lbar = 0.0;
        for (int i = 0; i "lt" N; ++i)
          Lbar += L[i] * weights[i];

        if (Lbar "gte" 0.5)  // bad learner
          continue;

        // 6. compute beta (low beta = high confidence)
        this.learners.Add(lrnr);  // lrnr is good
        double beta = Lbar / (1.0 - Lbar);
        this.betas.Add(beta); // each learner has a beta

        // 7. update training weights
        double sum = 0.0;
        for (int i = 0; i "lt" N; ++i)
        {
          weights[i] *= Math.Pow(beta, (1.0 - Lbar));
          sum += weights[i];
        }
        for (int i = 0; i "lt" N; ++i)
          weights[i] /= sum;

      } // t

    } // Train()

    public double Predict(double[] x)
    {
      // each learner makes a prediction
      double[] preds = new double[this.learners.Count];
      for (int i = 0; i "lt" this.learners.Count; ++i)
      {
        // preds[i] = lrnRate * this.learners[i].Predict(x);
        // scikit uses a lrnRate param, default val = 1.0
        preds[i] = this.learners[i].Predict(x);
      }

      // each learner has a "model weight" derived
      // from its beta-confidence weight
      double[] modelWts = new double[this.learners.Count];
      for (int i = 0; i "lt" this.learners.Count; ++i)
        modelWts[i] = Math.Log(1.0 / this.betas[i]);

      // final prediction is weighted median of 
      // the predictions
      double result = WeightedMedian(preds, modelWts);
      //double result = WeightedMean(preds, modelWts);

      return result;
    } // Predict()

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

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

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

    // ------------------------------------------------------
    // helper functions for AdaBoostRegressor
    // ------------------------------------------------------

    //private static double WeightedMean(double[] values,
    //  double[] weights)
    //{
    //  int n = values.Length;
    //  double sumWts = 0.0;
    //  double sum = 0.0;
    //  for (int i = 0; i "lt" n; ++i)
    //  {
    //    sumWts += weights[i];
    //    sum += values[i] * weights[i];
    //  }
    //  return sum / sumWts;
    //}

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

    private static double WeightedMedian(double[] values,
      double[] weights)
    {
      // no interpolation for even n
      // don't assume weights sum to 1.0
      int n = values.Length;
      double sumWts = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sumWts += weights[i];
      double thresh = sumWts / 2;
      int[] sortedIdxs = ArgSort(values);

      double accum = 0.0;
      for (int j = 0; j "lt" n; ++j)
      {
        accum += weights[sortedIdxs[j]];
        if (accum "gte" thresh)
          return values[sortedIdxs[j]];
      }
      return values[sortedIdxs[n - 1]];
    }

    // helper for WeightedMedian()
    private static int[] ArgSort(double[] values)
    {
      int n = values.Length;
      double[] copy = new double[n];
      int[] indices = new int[n];
      for (int i = 0; i "lt" n; ++i)
      {
        copy[i] = values[i];
        indices[i] = i;
      }
      Array.Sort(copy, indices);  // in parallel
      return indices;
    }

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

    private int SelectIdx(double[] probs)
    {
      // roulette wheel selection
      double p = this.rnd.NextDouble();
      double accum = 0.0;
      for (int i = 0; i "lt" probs.Length; ++i)
      {
        accum += probs[i];
        if (p "lt" accum)
          return i;
      }
      return probs.Length - 1;
    }

    private int[] SelectIdxs(double[] probs, int n)
    {
      // with replacement
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = this.SelectIdx(probs);
      return result;
    }

    private static double[][] ExtractRows(double[][] mat,
      int[] idxs)
    {
      int nRows = idxs.Length;  // of the result
      int nCols = mat[0].Length;  // source and result
      double[][] result = MatMake(nRows, nCols);

      for (int i = 0; i "lt" nRows; ++i)
      {
        for (int j = 0; j "lt" nCols; ++j)
        {
          int srcRow = idxs[i];
          result[i][j] = mat[srcRow][j];
        }
      }
      return result;
    }

    private static double[] ExtractVals(double[] vec,
      int[] idxs)
    {
      int nVals = idxs.Length;
      double[] result = new double[nVals];
      for (int i = 0; i "lt" nVals; ++i)
      {
        int srcIdx = idxs[i];
        result[i] = vec[srcIdx];
      }
      return result;
    }

  } // class AdaBoostRegressor

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

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

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

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

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

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

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

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

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

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

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

      // optionally delete rows in each Node to save space
      // when tree is part of an ensemble
      for (int i = 0; i "lt" this.tree.Count; ++i)
        if (this.tree[i] != null)
          this.tree[i].rows = null;
    }

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

    public double Predict(double[] x)
    {
      int p = 0;
      Node currNode = this.tree[p];
      while (currNode != null &&
        currNode.isLeaf == false &&
        p "lt" this.tree.Count)
      {
        if (x[currNode.colIdx] "lte" currNode.thresh)
          p = currNode.left;
        else
          p = currNode.right;
        currNode = this.tree[p];
      }
      return this.tree[p].value;
    }

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

    // helpers: MakeTree(), BestSplit(),
    // TreeTargetMean(), TreeTargetVariance()

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

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

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

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

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

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

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

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

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

        // construct the children, except colIdx and thresh
        Node leftNode = new Node();
        Node rightNode = new Node();

        // construct the children rows using split info
        // all info except colIdx and thresh
        List"lt"int"gt" leftIdxs = new List"lt"int"gt"();
        List"lt"int"gt" rightIdxs = new List"lt"int"gt"();
        for (int k = 0; k "lt" currNode.rows.Count; ++k)
        {
          int r = currNode.rows[k];
          if (this.trainX[r][colIdx] "lte" splitVal)
            leftIdxs.Add(r);
          else
            rightIdxs.Add(r);
        }

        // assign -1 to children if out of range
        leftNode.id = currNode.id * 2 + 1;
        if (leftNode.id "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.id = -1;
        leftNode.left = leftNode.id * 2 + 1;
        if (leftNode.left "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.left = -1;
        leftNode.right = leftNode.id * 2 + 2;
        if (leftNode.right "gt" (int)Math.Pow(2,
          (maxDepth + 1)) - 2) leftNode.right = -1;
        leftNode.rows = leftIdxs;
        leftNode.value =
          this.TreeTargetMean(leftNode.rows);
        this.tree[leftNode.id] = leftNode;

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

      } // i
      return;
    }

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

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

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

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

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

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

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

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

          // get row idxs where x is lte, gt thresh
          List"lt"int"gt" leftIdxs = new List"lt"int"gt"();
          List"lt"int"gt" rightIdxs = new List"lt"int"gt"();
          for (int k = 0; k "lt" nRows; ++k)
          {
            if (this.trainX[rows[k]][colIdx] "lte" thresh)
              leftIdxs.Add(rows[k]);
            else
              rightIdxs.Add(rows[k]);
          }

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

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

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

        } // each row
      } // j each col

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

    } // BestSplit()

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

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

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

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

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

  } // class DecisionTreeRegressor

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

} // ns

Training data:


# diabetes_norm_train_342.txt
# cols [0] to [9] predictors. col [10] target
# norm division constants:
# 100, -1, 100, 1000, 1000, 1000, 100, 10, 10, 1000, 1000
#
0.5900, 1.0000, 0.3210, 0.1010, 0.1570, 0.0932, 0.3800, 0.4000, 0.4860, 0.0870, 0.1510
0.4800, 0.0000, 0.2160, 0.0870, 0.1830, 0.1032, 0.7000, 0.3000, 0.3892, 0.0690, 0.0750
0.7200, 1.0000, 0.3050, 0.0930, 0.1560, 0.0936, 0.4100, 0.4000, 0.4673, 0.0850, 0.1410
0.2400, 0.0000, 0.2530, 0.0840, 0.1980, 0.1314, 0.4000, 0.5000, 0.4890, 0.0890, 0.2060
0.5000, 0.0000, 0.2300, 0.1010, 0.1920, 0.1254, 0.5200, 0.4000, 0.4291, 0.0800, 0.1350
0.2300, 0.0000, 0.2260, 0.0890, 0.1390, 0.0648, 0.6100, 0.2000, 0.4190, 0.0680, 0.0970
0.3600, 1.0000, 0.2200, 0.0900, 0.1600, 0.0996, 0.5000, 0.3000, 0.3951, 0.0820, 0.1380
0.6600, 1.0000, 0.2620, 0.1140, 0.2550, 0.1850, 0.5600, 0.4550, 0.4249, 0.0920, 0.0630
0.6000, 1.0000, 0.3210, 0.0830, 0.1790, 0.1194, 0.4200, 0.4000, 0.4477, 0.0940, 0.1100
0.2900, 0.0000, 0.3000, 0.0850, 0.1800, 0.0934, 0.4300, 0.4000, 0.5385, 0.0880, 0.3100
0.2200, 0.0000, 0.1860, 0.0970, 0.1140, 0.0576, 0.4600, 0.2000, 0.3951, 0.0830, 0.1010
0.5600, 1.0000, 0.2800, 0.0850, 0.1840, 0.1448, 0.3200, 0.6000, 0.3584, 0.0770, 0.0690
0.5300, 0.0000, 0.2370, 0.0920, 0.1860, 0.1092, 0.6200, 0.3000, 0.4304, 0.0810, 0.1790
0.5000, 1.0000, 0.2620, 0.0970, 0.1860, 0.1054, 0.4900, 0.4000, 0.5063, 0.0880, 0.1850
0.6100, 0.0000, 0.2400, 0.0910, 0.2020, 0.1154, 0.7200, 0.3000, 0.4291, 0.0730, 0.1180
0.3400, 1.0000, 0.2470, 0.1180, 0.2540, 0.1842, 0.3900, 0.7000, 0.5037, 0.0810, 0.1710
0.4700, 0.0000, 0.3030, 0.1090, 0.2070, 0.1002, 0.7000, 0.3000, 0.5215, 0.0980, 0.1660
0.6800, 1.0000, 0.2750, 0.1110, 0.2140, 0.1470, 0.3900, 0.5000, 0.4942, 0.0910, 0.1440
0.3800, 0.0000, 0.2540, 0.0840, 0.1620, 0.1030, 0.4200, 0.4000, 0.4443, 0.0870, 0.0970
0.4100, 0.0000, 0.2470, 0.0830, 0.1870, 0.1082, 0.6000, 0.3000, 0.4543, 0.0780, 0.1680
0.3500, 0.0000, 0.2110, 0.0820, 0.1560, 0.0878, 0.5000, 0.3000, 0.4511, 0.0950, 0.0680
0.2500, 1.0000, 0.2430, 0.0950, 0.1620, 0.0986, 0.5400, 0.3000, 0.3850, 0.0870, 0.0490
0.2500, 0.0000, 0.2600, 0.0920, 0.1870, 0.1204, 0.5600, 0.3000, 0.3970, 0.0880, 0.0680
0.6100, 1.0000, 0.3200, 0.1037, 0.2100, 0.0852, 0.3500, 0.6000, 0.6107, 0.1240, 0.2450
0.3100, 0.0000, 0.2970, 0.0880, 0.1670, 0.1034, 0.4800, 0.4000, 0.4357, 0.0780, 0.1840
0.3000, 1.0000, 0.2520, 0.0830, 0.1780, 0.1184, 0.3400, 0.5000, 0.4852, 0.0830, 0.2020
0.1900, 0.0000, 0.1920, 0.0870, 0.1240, 0.0540, 0.5700, 0.2000, 0.4174, 0.0900, 0.1370
0.4200, 0.0000, 0.3190, 0.0830, 0.1580, 0.0876, 0.5300, 0.3000, 0.4466, 0.1010, 0.0850
0.6300, 0.0000, 0.2440, 0.0730, 0.1600, 0.0914, 0.4800, 0.3000, 0.4635, 0.0780, 0.1310
0.6700, 1.0000, 0.2580, 0.1130, 0.1580, 0.0542, 0.6400, 0.2000, 0.5293, 0.1040, 0.2830
0.3200, 0.0000, 0.3050, 0.0890, 0.1820, 0.1106, 0.5600, 0.3000, 0.4344, 0.0890, 0.1290
0.4200, 0.0000, 0.2030, 0.0710, 0.1610, 0.0812, 0.6600, 0.2000, 0.4234, 0.0810, 0.0590
0.5800, 1.0000, 0.3800, 0.1030, 0.1500, 0.1072, 0.2200, 0.7000, 0.4644, 0.0980, 0.3410
0.5700, 0.0000, 0.2170, 0.0940, 0.1570, 0.0580, 0.8200, 0.2000, 0.4443, 0.0920, 0.0870
0.5300, 0.0000, 0.2050, 0.0780, 0.1470, 0.0842, 0.5200, 0.3000, 0.3989, 0.0750, 0.0650
0.6200, 1.0000, 0.2350, 0.0803, 0.2250, 0.1128, 0.8600, 0.2620, 0.4875, 0.0960, 0.1020
0.5200, 0.0000, 0.2850, 0.1100, 0.1950, 0.0972, 0.6000, 0.3000, 0.5242, 0.0850, 0.2650
0.4600, 0.0000, 0.2740, 0.0780, 0.1710, 0.0880, 0.5800, 0.3000, 0.4828, 0.0900, 0.2760
0.4800, 1.0000, 0.3300, 0.1230, 0.2530, 0.1636, 0.4400, 0.6000, 0.5425, 0.0970, 0.2520
0.4800, 1.0000, 0.2770, 0.0730, 0.1910, 0.1194, 0.4600, 0.4000, 0.4852, 0.0920, 0.0900
0.5000, 1.0000, 0.2560, 0.1010, 0.2290, 0.1622, 0.4300, 0.5000, 0.4779, 0.1140, 0.1000
0.2100, 0.0000, 0.2010, 0.0630, 0.1350, 0.0690, 0.5400, 0.3000, 0.4094, 0.0890, 0.0550
0.3200, 1.0000, 0.2540, 0.0903, 0.1530, 0.1004, 0.3400, 0.4500, 0.4533, 0.0830, 0.0610
0.5400, 0.0000, 0.2420, 0.0740, 0.2040, 0.1090, 0.8200, 0.2000, 0.4174, 0.1090, 0.0920
0.6100, 1.0000, 0.3270, 0.0970, 0.1770, 0.1184, 0.2900, 0.6000, 0.4997, 0.0870, 0.2590
0.5600, 1.0000, 0.2310, 0.1040, 0.1810, 0.1164, 0.4700, 0.4000, 0.4477, 0.0790, 0.0530
0.3300, 0.0000, 0.2530, 0.0850, 0.1550, 0.0850, 0.5100, 0.3000, 0.4554, 0.0700, 0.1900
0.2700, 0.0000, 0.1960, 0.0780, 0.1280, 0.0680, 0.4300, 0.3000, 0.4443, 0.0710, 0.1420
0.6700, 1.0000, 0.2250, 0.0980, 0.1910, 0.1192, 0.6100, 0.3000, 0.3989, 0.0860, 0.0750
0.3700, 1.0000, 0.2770, 0.0930, 0.1800, 0.1194, 0.3000, 0.6000, 0.5030, 0.0880, 0.1420
0.5800, 0.0000, 0.2570, 0.0990, 0.1570, 0.0916, 0.4900, 0.3000, 0.4407, 0.0930, 0.1550
0.6500, 1.0000, 0.2790, 0.1030, 0.1590, 0.0968, 0.4200, 0.4000, 0.4615, 0.0860, 0.2250
0.3400, 0.0000, 0.2550, 0.0930, 0.2180, 0.1440, 0.5700, 0.4000, 0.4443, 0.0880, 0.0590
0.4600, 0.0000, 0.2490, 0.1150, 0.1980, 0.1296, 0.5400, 0.4000, 0.4277, 0.1030, 0.1040
0.3500, 0.0000, 0.2870, 0.0970, 0.2040, 0.1268, 0.6400, 0.3000, 0.4190, 0.0930, 0.1820
0.3700, 0.0000, 0.2180, 0.0840, 0.1840, 0.1010, 0.7300, 0.3000, 0.3912, 0.0930, 0.1280
0.3700, 0.0000, 0.3020, 0.0870, 0.1660, 0.0960, 0.4000, 0.4150, 0.5011, 0.0870, 0.0520
0.4100, 0.0000, 0.2050, 0.0800, 0.1240, 0.0488, 0.6400, 0.2000, 0.4025, 0.0750, 0.0370
0.6000, 0.0000, 0.2040, 0.1050, 0.1980, 0.0784, 0.9900, 0.2000, 0.4635, 0.0790, 0.1700
0.6600, 1.0000, 0.2400, 0.0980, 0.2360, 0.1464, 0.5800, 0.4000, 0.5063, 0.0960, 0.1700
0.2900, 0.0000, 0.2600, 0.0830, 0.1410, 0.0652, 0.6400, 0.2000, 0.4078, 0.0830, 0.0610
0.3700, 1.0000, 0.2680, 0.0790, 0.1570, 0.0980, 0.2800, 0.6000, 0.5043, 0.0960, 0.1440
0.4100, 1.0000, 0.2570, 0.0830, 0.1810, 0.1066, 0.6600, 0.3000, 0.3738, 0.0850, 0.0520
0.3900, 0.0000, 0.2290, 0.0770, 0.2040, 0.1432, 0.4600, 0.4000, 0.4304, 0.0740, 0.1280
0.6700, 1.0000, 0.2400, 0.0830, 0.1430, 0.0772, 0.4900, 0.3000, 0.4431, 0.0940, 0.0710
0.3600, 1.0000, 0.2410, 0.1120, 0.1930, 0.1250, 0.3500, 0.6000, 0.5106, 0.0950, 0.1630
0.4600, 1.0000, 0.2470, 0.0850, 0.1740, 0.1232, 0.3000, 0.6000, 0.4644, 0.0960, 0.1500
0.6000, 1.0000, 0.2500, 0.0897, 0.1850, 0.1208, 0.4600, 0.4020, 0.4511, 0.0920, 0.0970
0.5900, 1.0000, 0.2360, 0.0830, 0.1650, 0.1000, 0.4700, 0.4000, 0.4500, 0.0920, 0.1600
0.5300, 0.0000, 0.2210, 0.0930, 0.1340, 0.0762, 0.4600, 0.3000, 0.4078, 0.0960, 0.1780
0.4800, 0.0000, 0.1990, 0.0910, 0.1890, 0.1096, 0.6900, 0.3000, 0.3951, 0.1010, 0.0480
0.4800, 0.0000, 0.2950, 0.1310, 0.2070, 0.1322, 0.4700, 0.4000, 0.4935, 0.1060, 0.2700
0.6600, 1.0000, 0.2600, 0.0910, 0.2640, 0.1466, 0.6500, 0.4000, 0.5568, 0.0870, 0.2020
0.5200, 1.0000, 0.2450, 0.0940, 0.2170, 0.1494, 0.4800, 0.5000, 0.4585, 0.0890, 0.1110
0.5200, 1.0000, 0.2660, 0.1110, 0.2090, 0.1264, 0.6100, 0.3000, 0.4682, 0.1090, 0.0850
0.4600, 1.0000, 0.2350, 0.0870, 0.1810, 0.1148, 0.4400, 0.4000, 0.4710, 0.0980, 0.0420
0.4000, 1.0000, 0.2900, 0.1150, 0.0970, 0.0472, 0.3500, 0.2770, 0.4304, 0.0950, 0.1700
0.2200, 0.0000, 0.2300, 0.0730, 0.1610, 0.0978, 0.5400, 0.3000, 0.3829, 0.0910, 0.2000
0.5000, 0.0000, 0.2100, 0.0880, 0.1400, 0.0718, 0.3500, 0.4000, 0.5112, 0.0710, 0.2520
0.2000, 0.0000, 0.2290, 0.0870, 0.1910, 0.1282, 0.5300, 0.4000, 0.3892, 0.0850, 0.1130
0.6800, 0.0000, 0.2750, 0.1070, 0.2410, 0.1496, 0.6400, 0.4000, 0.4920, 0.0900, 0.1430
0.5200, 1.0000, 0.2430, 0.0860, 0.1970, 0.1336, 0.4400, 0.5000, 0.4575, 0.0910, 0.0510
0.4400, 0.0000, 0.2310, 0.0870, 0.2130, 0.1264, 0.7700, 0.3000, 0.3871, 0.0720, 0.0520
0.3800, 0.0000, 0.2730, 0.0810, 0.1460, 0.0816, 0.4700, 0.3000, 0.4466, 0.0810, 0.2100
0.4900, 0.0000, 0.2270, 0.0653, 0.1680, 0.0962, 0.6200, 0.2710, 0.3892, 0.0600, 0.0650
0.6100, 0.0000, 0.3300, 0.0950, 0.1820, 0.1148, 0.5400, 0.3000, 0.4190, 0.0740, 0.1410
0.2900, 1.0000, 0.1940, 0.0830, 0.1520, 0.1058, 0.3900, 0.4000, 0.3584, 0.0830, 0.0550
0.6100, 0.0000, 0.2580, 0.0980, 0.2350, 0.1258, 0.7600, 0.3000, 0.5112, 0.0820, 0.1340
0.3400, 1.0000, 0.2260, 0.0750, 0.1660, 0.0918, 0.6000, 0.3000, 0.4263, 0.1080, 0.0420
0.3600, 0.0000, 0.2190, 0.0890, 0.1890, 0.1052, 0.6800, 0.3000, 0.4369, 0.0960, 0.1110
0.5200, 0.0000, 0.2400, 0.0830, 0.1670, 0.0866, 0.7100, 0.2000, 0.3850, 0.0940, 0.0980
0.6100, 0.0000, 0.3120, 0.0790, 0.2350, 0.1568, 0.4700, 0.5000, 0.5050, 0.0960, 0.1640
0.4300, 0.0000, 0.2680, 0.1230, 0.1930, 0.1022, 0.6700, 0.3000, 0.4779, 0.0940, 0.0480
0.3500, 0.0000, 0.2040, 0.0650, 0.1870, 0.1056, 0.6700, 0.2790, 0.4277, 0.0780, 0.0960
0.2700, 0.0000, 0.2480, 0.0910, 0.1890, 0.1068, 0.6900, 0.3000, 0.4190, 0.0690, 0.0900
0.2900, 0.0000, 0.2100, 0.0710, 0.1560, 0.0970, 0.3800, 0.4000, 0.4654, 0.0900, 0.1620
0.6400, 1.0000, 0.2730, 0.1090, 0.1860, 0.1076, 0.3800, 0.5000, 0.5308, 0.0990, 0.1500
0.4100, 0.0000, 0.3460, 0.0873, 0.2050, 0.1426, 0.4100, 0.5000, 0.4673, 0.1100, 0.2790
0.4900, 1.0000, 0.2590, 0.0910, 0.1780, 0.1066, 0.5200, 0.3000, 0.4575, 0.0750, 0.0920
0.4800, 0.0000, 0.2040, 0.0980, 0.2090, 0.1394, 0.4600, 0.5000, 0.4771, 0.0780, 0.0830
0.5300, 0.0000, 0.2800, 0.0880, 0.2330, 0.1438, 0.5800, 0.4000, 0.5050, 0.0910, 0.1280
0.5300, 1.0000, 0.2220, 0.1130, 0.1970, 0.1152, 0.6700, 0.3000, 0.4304, 0.1000, 0.1020
0.2300, 0.0000, 0.2900, 0.0900, 0.2160, 0.1314, 0.6500, 0.3000, 0.4585, 0.0910, 0.3020
0.6500, 1.0000, 0.3020, 0.0980, 0.2190, 0.1606, 0.4000, 0.5000, 0.4522, 0.0840, 0.1980
0.4100, 0.0000, 0.3240, 0.0940, 0.1710, 0.1044, 0.5600, 0.3000, 0.3970, 0.0760, 0.0950
0.5500, 1.0000, 0.2340, 0.0830, 0.1660, 0.1016, 0.4600, 0.4000, 0.4522, 0.0960, 0.0530
0.2200, 0.0000, 0.1930, 0.0820, 0.1560, 0.0932, 0.5200, 0.3000, 0.3989, 0.0710, 0.1340
0.5600, 0.0000, 0.3100, 0.0787, 0.1870, 0.1414, 0.3400, 0.5500, 0.4060, 0.0900, 0.1440
0.5400, 1.0000, 0.3060, 0.1033, 0.1440, 0.0798, 0.3000, 0.4800, 0.5142, 0.1010, 0.2320
0.5900, 1.0000, 0.2550, 0.0953, 0.1900, 0.1394, 0.3500, 0.5430, 0.4357, 0.1170, 0.0810
0.6000, 1.0000, 0.2340, 0.0880, 0.1530, 0.0898, 0.5800, 0.3000, 0.3258, 0.0950, 0.1040
0.5400, 0.0000, 0.2680, 0.0870, 0.2060, 0.1220, 0.6800, 0.3000, 0.4382, 0.0800, 0.0590
0.2500, 0.0000, 0.2830, 0.0870, 0.1930, 0.1280, 0.4900, 0.4000, 0.4382, 0.0920, 0.2460
0.5400, 1.0000, 0.2770, 0.1130, 0.2000, 0.1284, 0.3700, 0.5000, 0.5153, 0.1130, 0.2970
0.5500, 0.0000, 0.3660, 0.1130, 0.1990, 0.0944, 0.4300, 0.4630, 0.5730, 0.0970, 0.2580
0.4000, 1.0000, 0.2650, 0.0930, 0.2360, 0.1470, 0.3700, 0.7000, 0.5561, 0.0920, 0.2290
0.6200, 1.0000, 0.3180, 0.1150, 0.1990, 0.1286, 0.4400, 0.5000, 0.4883, 0.0980, 0.2750
0.6500, 0.0000, 0.2440, 0.1200, 0.2220, 0.1356, 0.3700, 0.6000, 0.5509, 0.1240, 0.2810
0.3300, 1.0000, 0.2540, 0.1020, 0.2060, 0.1410, 0.3900, 0.5000, 0.4868, 0.1050, 0.1790
0.5300, 0.0000, 0.2200, 0.0940, 0.1750, 0.0880, 0.5900, 0.3000, 0.4942, 0.0980, 0.2000
0.3500, 0.0000, 0.2680, 0.0980, 0.1620, 0.1036, 0.4500, 0.4000, 0.4205, 0.0860, 0.2000
0.6600, 0.0000, 0.2800, 0.1010, 0.1950, 0.1292, 0.4000, 0.5000, 0.4860, 0.0940, 0.1730
0.6200, 1.0000, 0.3390, 0.1010, 0.2210, 0.1564, 0.3500, 0.6000, 0.4997, 0.1030, 0.1800
0.5000, 1.0000, 0.2960, 0.0943, 0.3000, 0.2424, 0.3300, 0.9090, 0.4812, 0.1090, 0.0840
0.4700, 0.0000, 0.2860, 0.0970, 0.1640, 0.0906, 0.5600, 0.3000, 0.4466, 0.0880, 0.1210
0.4700, 1.0000, 0.2560, 0.0940, 0.1650, 0.0748, 0.4000, 0.4000, 0.5526, 0.0930, 0.1610
0.2400, 0.0000, 0.2070, 0.0870, 0.1490, 0.0806, 0.6100, 0.2000, 0.3611, 0.0780, 0.0990
0.5800, 1.0000, 0.2620, 0.0910, 0.2170, 0.1242, 0.7100, 0.3000, 0.4691, 0.0680, 0.1090
0.3400, 0.0000, 0.2060, 0.0870, 0.1850, 0.1122, 0.5800, 0.3000, 0.4304, 0.0740, 0.1150
0.5100, 0.0000, 0.2790, 0.0960, 0.1960, 0.1222, 0.4200, 0.5000, 0.5069, 0.1200, 0.2680
0.3100, 1.0000, 0.3530, 0.1250, 0.1870, 0.1124, 0.4800, 0.4000, 0.4890, 0.1090, 0.2740
0.2200, 0.0000, 0.1990, 0.0750, 0.1750, 0.1086, 0.5400, 0.3000, 0.4127, 0.0720, 0.1580
0.5300, 1.0000, 0.2440, 0.0920, 0.2140, 0.1460, 0.5000, 0.4000, 0.4500, 0.0970, 0.1070
0.3700, 1.0000, 0.2140, 0.0830, 0.1280, 0.0696, 0.4900, 0.3000, 0.3850, 0.0840, 0.0830
0.2800, 0.0000, 0.3040, 0.0850, 0.1980, 0.1156, 0.6700, 0.3000, 0.4344, 0.0800, 0.1030
0.4700, 0.0000, 0.3160, 0.0840, 0.1540, 0.0880, 0.3000, 0.5100, 0.5199, 0.1050, 0.2720
0.2300, 0.0000, 0.1880, 0.0780, 0.1450, 0.0720, 0.6300, 0.2000, 0.3912, 0.0860, 0.0850
0.5000, 0.0000, 0.3100, 0.1230, 0.1780, 0.1050, 0.4800, 0.4000, 0.4828, 0.0880, 0.2800
0.5800, 1.0000, 0.3670, 0.1170, 0.1660, 0.0938, 0.4400, 0.4000, 0.4949, 0.1090, 0.3360
0.5500, 0.0000, 0.3210, 0.1100, 0.1640, 0.0842, 0.4200, 0.4000, 0.5242, 0.0900, 0.2810
0.6000, 1.0000, 0.2770, 0.1070, 0.1670, 0.1146, 0.3800, 0.4000, 0.4277, 0.0950, 0.1180
0.4100, 0.0000, 0.3080, 0.0810, 0.2140, 0.1520, 0.2800, 0.7600, 0.5136, 0.1230, 0.3170
0.6000, 1.0000, 0.2750, 0.1060, 0.2290, 0.1438, 0.5100, 0.4000, 0.5142, 0.0910, 0.2350
0.4000, 0.0000, 0.2690, 0.0920, 0.2030, 0.1198, 0.7000, 0.3000, 0.4190, 0.0810, 0.0600
0.5700, 1.0000, 0.3070, 0.0900, 0.2040, 0.1478, 0.3400, 0.6000, 0.4710, 0.0930, 0.1740
0.3700, 0.0000, 0.3830, 0.1130, 0.1650, 0.0946, 0.5300, 0.3000, 0.4466, 0.0790, 0.2590
0.4000, 1.0000, 0.3190, 0.0950, 0.1980, 0.1356, 0.3800, 0.5000, 0.4804, 0.0930, 0.1780
0.3300, 0.0000, 0.3500, 0.0890, 0.2000, 0.1304, 0.4200, 0.4760, 0.4927, 0.1010, 0.1280
0.3200, 1.0000, 0.2780, 0.0890, 0.2160, 0.1462, 0.5500, 0.4000, 0.4304, 0.0910, 0.0960
0.3500, 1.0000, 0.2590, 0.0810, 0.1740, 0.1024, 0.3100, 0.6000, 0.5313, 0.0820, 0.1260
0.5500, 0.0000, 0.3290, 0.1020, 0.1640, 0.1062, 0.4100, 0.4000, 0.4431, 0.0890, 0.2880
0.4900, 0.0000, 0.2600, 0.0930, 0.1830, 0.1002, 0.6400, 0.3000, 0.4543, 0.0880, 0.0880
0.3900, 1.0000, 0.2630, 0.1150, 0.2180, 0.1582, 0.3200, 0.7000, 0.4935, 0.1090, 0.2920
0.6000, 1.0000, 0.2230, 0.1130, 0.1860, 0.1258, 0.4600, 0.4000, 0.4263, 0.0940, 0.0710
0.6700, 1.0000, 0.2830, 0.0930, 0.2040, 0.1322, 0.4900, 0.4000, 0.4736, 0.0920, 0.1970
0.4100, 1.0000, 0.3200, 0.1090, 0.2510, 0.1706, 0.4900, 0.5000, 0.5056, 0.1030, 0.1860
0.4400, 0.0000, 0.2540, 0.0950, 0.1620, 0.0926, 0.5300, 0.3000, 0.4407, 0.0830, 0.0250
0.4800, 1.0000, 0.2330, 0.0893, 0.2120, 0.1428, 0.4600, 0.4610, 0.4754, 0.0980, 0.0840
0.4500, 0.0000, 0.2030, 0.0743, 0.1900, 0.1262, 0.4900, 0.3880, 0.4304, 0.0790, 0.0960
0.4700, 0.0000, 0.3040, 0.1200, 0.1990, 0.1200, 0.4600, 0.4000, 0.5106, 0.0870, 0.1950
0.4600, 0.0000, 0.2060, 0.0730, 0.1720, 0.1070, 0.5100, 0.3000, 0.4249, 0.0800, 0.0530
0.3600, 1.0000, 0.3230, 0.1150, 0.2860, 0.1994, 0.3900, 0.7000, 0.5472, 0.1120, 0.2170
0.3400, 0.0000, 0.2920, 0.0730, 0.1720, 0.1082, 0.4900, 0.4000, 0.4304, 0.0910, 0.1720
0.5300, 1.0000, 0.3310, 0.1170, 0.1830, 0.1190, 0.4800, 0.4000, 0.4382, 0.1060, 0.1310
0.6100, 0.0000, 0.2460, 0.1010, 0.2090, 0.1068, 0.7700, 0.3000, 0.4836, 0.0880, 0.2140
0.3700, 0.0000, 0.2020, 0.0810, 0.1620, 0.0878, 0.6300, 0.3000, 0.4025, 0.0880, 0.0590
0.3300, 1.0000, 0.2080, 0.0840, 0.1250, 0.0702, 0.4600, 0.3000, 0.3784, 0.0660, 0.0700
0.6800, 0.0000, 0.3280, 0.1057, 0.2050, 0.1164, 0.4000, 0.5130, 0.5493, 0.1170, 0.2200
0.4900, 1.0000, 0.3190, 0.0940, 0.2340, 0.1558, 0.3400, 0.7000, 0.5398, 0.1220, 0.2680
0.4800, 0.0000, 0.2390, 0.1090, 0.2320, 0.1052, 0.3700, 0.6000, 0.6107, 0.0960, 0.1520
0.5500, 1.0000, 0.2450, 0.0840, 0.1790, 0.1058, 0.6600, 0.3000, 0.3584, 0.0870, 0.0470
0.4300, 0.0000, 0.2210, 0.0660, 0.1340, 0.0772, 0.4500, 0.3000, 0.4078, 0.0800, 0.0740
0.6000, 1.0000, 0.3300, 0.0970, 0.2170, 0.1256, 0.4500, 0.5000, 0.5447, 0.1120, 0.2950
0.3100, 1.0000, 0.1900, 0.0930, 0.1370, 0.0730, 0.4700, 0.3000, 0.4443, 0.0780, 0.1010
0.5300, 1.0000, 0.2730, 0.0820, 0.1190, 0.0550, 0.3900, 0.3000, 0.4828, 0.0930, 0.1510
0.6700, 0.0000, 0.2280, 0.0870, 0.1660, 0.0986, 0.5200, 0.3000, 0.4344, 0.0920, 0.1270
0.6100, 1.0000, 0.2820, 0.1060, 0.2040, 0.1320, 0.5200, 0.4000, 0.4605, 0.0960, 0.2370
0.6200, 0.0000, 0.2890, 0.0873, 0.2060, 0.1272, 0.3300, 0.6240, 0.5434, 0.0990, 0.2250
0.6000, 0.0000, 0.2560, 0.0870, 0.2070, 0.1258, 0.6900, 0.3000, 0.4111, 0.0840, 0.0810
0.4200, 0.0000, 0.2490, 0.0910, 0.2040, 0.1418, 0.3800, 0.5000, 0.4796, 0.0890, 0.1510
0.3800, 1.0000, 0.2680, 0.1050, 0.1810, 0.1192, 0.3700, 0.5000, 0.4820, 0.0910, 0.1070
0.6200, 0.0000, 0.2240, 0.0790, 0.2220, 0.1474, 0.5900, 0.4000, 0.4357, 0.0760, 0.0640
0.6100, 1.0000, 0.2690, 0.1110, 0.2360, 0.1724, 0.3900, 0.6000, 0.4812, 0.0890, 0.1380
0.6100, 1.0000, 0.2310, 0.1130, 0.1860, 0.1144, 0.4700, 0.4000, 0.4812, 0.1050, 0.1850
0.5300, 0.0000, 0.2860, 0.0880, 0.1710, 0.0988, 0.4100, 0.4000, 0.5050, 0.0990, 0.2650
0.2800, 1.0000, 0.2470, 0.0970, 0.1750, 0.0996, 0.3200, 0.5000, 0.5380, 0.0870, 0.1010
0.2600, 1.0000, 0.3030, 0.0890, 0.2180, 0.1522, 0.3100, 0.7000, 0.5159, 0.0820, 0.1370
0.3000, 0.0000, 0.2130, 0.0870, 0.1340, 0.0630, 0.6300, 0.2000, 0.3689, 0.0660, 0.1430
0.5000, 0.0000, 0.2610, 0.1090, 0.2430, 0.1606, 0.6200, 0.4000, 0.4625, 0.0890, 0.1410
0.4800, 0.0000, 0.2020, 0.0950, 0.1870, 0.1174, 0.5300, 0.4000, 0.4419, 0.0850, 0.0790
0.5100, 0.0000, 0.2520, 0.1030, 0.1760, 0.1122, 0.3700, 0.5000, 0.4898, 0.0900, 0.2920
0.4700, 1.0000, 0.2250, 0.0820, 0.1310, 0.0668, 0.4100, 0.3000, 0.4754, 0.0890, 0.1780
0.6400, 1.0000, 0.2350, 0.0970, 0.2030, 0.1290, 0.5900, 0.3000, 0.4318, 0.0770, 0.0910
0.5100, 1.0000, 0.2590, 0.0760, 0.2400, 0.1690, 0.3900, 0.6000, 0.5075, 0.0960, 0.1160
0.3000, 0.0000, 0.2090, 0.1040, 0.1520, 0.0838, 0.4700, 0.3000, 0.4663, 0.0970, 0.0860
0.5600, 1.0000, 0.2870, 0.0990, 0.2080, 0.1464, 0.3900, 0.5000, 0.4727, 0.0970, 0.1220
0.4200, 0.0000, 0.2210, 0.0850, 0.2130, 0.1386, 0.6000, 0.4000, 0.4277, 0.0940, 0.0720
0.6200, 1.0000, 0.2670, 0.1150, 0.1830, 0.1240, 0.3500, 0.5000, 0.4788, 0.1000, 0.1290
0.3400, 0.0000, 0.3140, 0.0870, 0.1490, 0.0938, 0.4600, 0.3000, 0.3829, 0.0770, 0.1420
0.6000, 0.0000, 0.2220, 0.1047, 0.2210, 0.1054, 0.6000, 0.3680, 0.5628, 0.0930, 0.0900
0.6400, 0.0000, 0.2100, 0.0923, 0.2270, 0.1468, 0.6500, 0.3490, 0.4331, 0.1020, 0.1580
0.3900, 1.0000, 0.2120, 0.0900, 0.1820, 0.1104, 0.6000, 0.3000, 0.4060, 0.0980, 0.0390
0.7100, 1.0000, 0.2650, 0.1050, 0.2810, 0.1736, 0.5500, 0.5000, 0.5568, 0.0840, 0.1960
0.4800, 1.0000, 0.2920, 0.1100, 0.2180, 0.1516, 0.3900, 0.6000, 0.4920, 0.0980, 0.2220
0.7900, 1.0000, 0.2700, 0.1030, 0.1690, 0.1108, 0.3700, 0.5000, 0.4663, 0.1100, 0.2770
0.4000, 0.0000, 0.3070, 0.0990, 0.1770, 0.0854, 0.5000, 0.4000, 0.5338, 0.0850, 0.0990
0.4900, 1.0000, 0.2880, 0.0920, 0.2070, 0.1400, 0.4400, 0.5000, 0.4745, 0.0920, 0.1960
0.5100, 0.0000, 0.3060, 0.1030, 0.1980, 0.1066, 0.5700, 0.3000, 0.5148, 0.1000, 0.2020
0.5700, 0.0000, 0.3010, 0.1170, 0.2020, 0.1396, 0.4200, 0.5000, 0.4625, 0.1200, 0.1550
0.5900, 1.0000, 0.2470, 0.1140, 0.1520, 0.1048, 0.2900, 0.5000, 0.4511, 0.0880, 0.0770
0.5100, 0.0000, 0.2770, 0.0990, 0.2290, 0.1456, 0.6900, 0.3000, 0.4277, 0.0770, 0.1910
0.7400, 0.0000, 0.2980, 0.1010, 0.1710, 0.1048, 0.5000, 0.3000, 0.4394, 0.0860, 0.0700
0.6700, 0.0000, 0.2670, 0.1050, 0.2250, 0.1354, 0.6900, 0.3000, 0.4635, 0.0960, 0.0730
0.4900, 0.0000, 0.1980, 0.0880, 0.1880, 0.1148, 0.5700, 0.3000, 0.4394, 0.0930, 0.0490
0.5700, 0.0000, 0.2330, 0.0880, 0.1550, 0.0636, 0.7800, 0.2000, 0.4205, 0.0780, 0.0650
0.5600, 1.0000, 0.3510, 0.1230, 0.1640, 0.0950, 0.3800, 0.4000, 0.5043, 0.1170, 0.2630
0.5200, 1.0000, 0.2970, 0.1090, 0.2280, 0.1628, 0.3100, 0.8000, 0.5142, 0.1030, 0.2480
0.6900, 0.0000, 0.2930, 0.1240, 0.2230, 0.1390, 0.5400, 0.4000, 0.5011, 0.1020, 0.2960
0.3700, 0.0000, 0.2030, 0.0830, 0.1850, 0.1246, 0.3800, 0.5000, 0.4719, 0.0880, 0.2140
0.2400, 0.0000, 0.2250, 0.0890, 0.1410, 0.0680, 0.5200, 0.3000, 0.4654, 0.0840, 0.1850
0.5500, 1.0000, 0.2270, 0.0930, 0.1540, 0.0942, 0.5300, 0.3000, 0.3526, 0.0750, 0.0780
0.3600, 0.0000, 0.2280, 0.0870, 0.1780, 0.1160, 0.4100, 0.4000, 0.4654, 0.0820, 0.0930
0.4200, 1.0000, 0.2400, 0.1070, 0.1500, 0.0850, 0.4400, 0.3000, 0.4654, 0.0960, 0.2520
0.2100, 0.0000, 0.2420, 0.0760, 0.1470, 0.0770, 0.5300, 0.3000, 0.4443, 0.0790, 0.1500
0.4100, 0.0000, 0.2020, 0.0620, 0.1530, 0.0890, 0.5000, 0.3000, 0.4249, 0.0890, 0.0770
0.5700, 1.0000, 0.2940, 0.1090, 0.1600, 0.0876, 0.3100, 0.5000, 0.5333, 0.0920, 0.2080
0.2000, 1.0000, 0.2210, 0.0870, 0.1710, 0.0996, 0.5800, 0.3000, 0.4205, 0.0780, 0.0770
0.6700, 1.0000, 0.2360, 0.1113, 0.1890, 0.1054, 0.7000, 0.2700, 0.4220, 0.0930, 0.1080
0.3400, 0.0000, 0.2520, 0.0770, 0.1890, 0.1206, 0.5300, 0.4000, 0.4344, 0.0790, 0.1600
0.4100, 1.0000, 0.2490, 0.0860, 0.1920, 0.1150, 0.6100, 0.3000, 0.4382, 0.0940, 0.0530
0.3800, 1.0000, 0.3300, 0.0780, 0.3010, 0.2150, 0.5000, 0.6020, 0.5193, 0.1080, 0.2200
0.5100, 0.0000, 0.2350, 0.1010, 0.1950, 0.1210, 0.5100, 0.4000, 0.4745, 0.0940, 0.1540
0.5200, 1.0000, 0.2640, 0.0913, 0.2180, 0.1520, 0.3900, 0.5590, 0.4905, 0.0990, 0.2590
0.6700, 0.0000, 0.2980, 0.0800, 0.1720, 0.0934, 0.6300, 0.3000, 0.4357, 0.0820, 0.0900
0.6100, 0.0000, 0.3000, 0.1080, 0.1940, 0.1000, 0.5200, 0.3730, 0.5347, 0.1050, 0.2460
0.6700, 1.0000, 0.2500, 0.1117, 0.1460, 0.0934, 0.3300, 0.4420, 0.4585, 0.1030, 0.1240
0.5600, 0.0000, 0.2700, 0.1050, 0.2470, 0.1606, 0.5400, 0.5000, 0.5088, 0.0940, 0.0670
0.6400, 0.0000, 0.2000, 0.0747, 0.1890, 0.1148, 0.6200, 0.3050, 0.4111, 0.0910, 0.0720
0.5800, 1.0000, 0.2550, 0.1120, 0.1630, 0.1106, 0.2900, 0.6000, 0.4762, 0.0860, 0.2570
0.5500, 0.0000, 0.2820, 0.0910, 0.2500, 0.1402, 0.6700, 0.4000, 0.5366, 0.1030, 0.2620
0.6200, 1.0000, 0.3330, 0.1140, 0.1820, 0.1140, 0.3800, 0.5000, 0.5011, 0.0960, 0.2750
0.5700, 1.0000, 0.2560, 0.0960, 0.2000, 0.1330, 0.5200, 0.3850, 0.4318, 0.1050, 0.1770
0.2000, 1.0000, 0.2420, 0.0880, 0.1260, 0.0722, 0.4500, 0.3000, 0.3784, 0.0740, 0.0710
0.5300, 1.0000, 0.2210, 0.0980, 0.1650, 0.1052, 0.4700, 0.4000, 0.4159, 0.0810, 0.0470
0.3200, 1.0000, 0.3140, 0.0890, 0.1530, 0.0842, 0.5600, 0.3000, 0.4159, 0.0900, 0.1870
0.4100, 0.0000, 0.2310, 0.0860, 0.1480, 0.0780, 0.5800, 0.3000, 0.4094, 0.0600, 0.1250
0.6000, 0.0000, 0.2340, 0.0767, 0.2470, 0.1480, 0.6500, 0.3800, 0.5136, 0.0770, 0.0780
0.2600, 0.0000, 0.1880, 0.0830, 0.1910, 0.1036, 0.6900, 0.3000, 0.4522, 0.0690, 0.0510
0.3700, 0.0000, 0.3080, 0.1120, 0.2820, 0.1972, 0.4300, 0.7000, 0.5342, 0.1010, 0.2580
0.4500, 0.0000, 0.3200, 0.1100, 0.2240, 0.1342, 0.4500, 0.5000, 0.5412, 0.0930, 0.2150
0.6700, 0.0000, 0.3160, 0.1160, 0.1790, 0.0904, 0.4100, 0.4000, 0.5472, 0.1000, 0.3030
0.3400, 1.0000, 0.3550, 0.1200, 0.2330, 0.1466, 0.3400, 0.7000, 0.5568, 0.1010, 0.2430
0.5000, 0.0000, 0.3190, 0.0783, 0.2070, 0.1492, 0.3800, 0.5450, 0.4595, 0.0840, 0.0910
0.7100, 0.0000, 0.2950, 0.0970, 0.2270, 0.1516, 0.4500, 0.5000, 0.5024, 0.1080, 0.1500
0.5700, 1.0000, 0.3160, 0.1170, 0.2250, 0.1076, 0.4000, 0.6000, 0.5958, 0.1130, 0.3100
0.4900, 0.0000, 0.2030, 0.0930, 0.1840, 0.1030, 0.6100, 0.3000, 0.4605, 0.0930, 0.1530
0.3500, 0.0000, 0.4130, 0.0810, 0.1680, 0.1028, 0.3700, 0.5000, 0.4949, 0.0940, 0.3460
0.4100, 1.0000, 0.2120, 0.1020, 0.1840, 0.1004, 0.6400, 0.3000, 0.4585, 0.0790, 0.0630
0.7000, 1.0000, 0.2410, 0.0823, 0.1940, 0.1492, 0.3100, 0.6260, 0.4234, 0.1050, 0.0890
0.5200, 0.0000, 0.2300, 0.1070, 0.1790, 0.1237, 0.4250, 0.4210, 0.4159, 0.0930, 0.0500
0.6000, 0.0000, 0.2560, 0.0780, 0.1950, 0.0954, 0.9100, 0.2000, 0.3761, 0.0870, 0.0390
0.6200, 0.0000, 0.2250, 0.1250, 0.2150, 0.0990, 0.9800, 0.2000, 0.4500, 0.0950, 0.1030
0.4400, 1.0000, 0.3820, 0.1230, 0.2010, 0.1266, 0.4400, 0.5000, 0.5024, 0.0920, 0.3080
0.2800, 1.0000, 0.1920, 0.0810, 0.1550, 0.0946, 0.5100, 0.3000, 0.3850, 0.0870, 0.1160
0.5800, 1.0000, 0.2900, 0.0850, 0.1560, 0.1092, 0.3600, 0.4000, 0.3989, 0.0860, 0.1450
0.3900, 1.0000, 0.2400, 0.0897, 0.1900, 0.1136, 0.5200, 0.3650, 0.4804, 0.1010, 0.0740
0.3400, 1.0000, 0.2060, 0.0980, 0.1830, 0.0920, 0.8300, 0.2000, 0.3689, 0.0920, 0.0450
0.6500, 0.0000, 0.2630, 0.0700, 0.2440, 0.1662, 0.5100, 0.5000, 0.4898, 0.0980, 0.1150
0.6600, 1.0000, 0.3460, 0.1150, 0.2040, 0.1394, 0.3600, 0.6000, 0.4963, 0.1090, 0.2640
0.5100, 0.0000, 0.2340, 0.0870, 0.2200, 0.1088, 0.9300, 0.2000, 0.4511, 0.0820, 0.0870
0.5000, 1.0000, 0.2920, 0.1190, 0.1620, 0.0852, 0.5400, 0.3000, 0.4736, 0.0950, 0.2020
0.5900, 1.0000, 0.2720, 0.1070, 0.1580, 0.1020, 0.3900, 0.4000, 0.4443, 0.0930, 0.1270
0.5200, 0.0000, 0.2700, 0.0783, 0.1340, 0.0730, 0.4400, 0.3050, 0.4443, 0.0690, 0.1820
0.6900, 1.0000, 0.2450, 0.1080, 0.2430, 0.1364, 0.4000, 0.6000, 0.5808, 0.1000, 0.2410
0.5300, 0.0000, 0.2410, 0.1050, 0.1840, 0.1134, 0.4600, 0.4000, 0.4812, 0.0950, 0.0660
0.4700, 1.0000, 0.2530, 0.0980, 0.1730, 0.1056, 0.4400, 0.4000, 0.4762, 0.1080, 0.0940
0.5200, 0.0000, 0.2880, 0.1130, 0.2800, 0.1740, 0.6700, 0.4000, 0.5273, 0.0860, 0.2830
0.3900, 0.0000, 0.2090, 0.0950, 0.1500, 0.0656, 0.6800, 0.2000, 0.4407, 0.0950, 0.0640
0.6700, 1.0000, 0.2300, 0.0700, 0.1840, 0.1280, 0.3500, 0.5000, 0.4654, 0.0990, 0.1020
0.5900, 1.0000, 0.2410, 0.0960, 0.1700, 0.0986, 0.5400, 0.3000, 0.4466, 0.0850, 0.2000
0.5100, 1.0000, 0.2810, 0.1060, 0.2020, 0.1222, 0.5500, 0.4000, 0.4820, 0.0870, 0.2650
0.2300, 1.0000, 0.1800, 0.0780, 0.1710, 0.0960, 0.4800, 0.4000, 0.4905, 0.0920, 0.0940
0.6800, 0.0000, 0.2590, 0.0930, 0.2530, 0.1812, 0.5300, 0.5000, 0.4543, 0.0980, 0.2300
0.4400, 0.0000, 0.2150, 0.0850, 0.1570, 0.0922, 0.5500, 0.3000, 0.3892, 0.0840, 0.1810
0.6000, 1.0000, 0.2430, 0.1030, 0.1410, 0.0866, 0.3300, 0.4000, 0.4673, 0.0780, 0.1560
0.5200, 0.0000, 0.2450, 0.0900, 0.1980, 0.1290, 0.2900, 0.7000, 0.5298, 0.0860, 0.2330
0.3800, 0.0000, 0.2130, 0.0720, 0.1650, 0.0602, 0.8800, 0.2000, 0.4431, 0.0900, 0.0600
0.6100, 0.0000, 0.2580, 0.0900, 0.2800, 0.1954, 0.5500, 0.5000, 0.4997, 0.0900, 0.2190
0.6800, 1.0000, 0.2480, 0.1010, 0.2210, 0.1514, 0.6000, 0.4000, 0.3871, 0.0870, 0.0800
0.2800, 1.0000, 0.3150, 0.0830, 0.2280, 0.1494, 0.3800, 0.6000, 0.5313, 0.0830, 0.0680
0.6500, 1.0000, 0.3350, 0.1020, 0.1900, 0.1262, 0.3500, 0.5000, 0.4970, 0.1020, 0.3320
0.6900, 0.0000, 0.2810, 0.1130, 0.2340, 0.1428, 0.5200, 0.4000, 0.5278, 0.0770, 0.2480
0.5100, 0.0000, 0.2430, 0.0853, 0.1530, 0.0716, 0.7100, 0.2150, 0.3951, 0.0820, 0.0840
0.2900, 0.0000, 0.3500, 0.0983, 0.2040, 0.1426, 0.5000, 0.4080, 0.4043, 0.0910, 0.2000
0.5500, 1.0000, 0.2350, 0.0930, 0.1770, 0.1268, 0.4100, 0.4000, 0.3829, 0.0830, 0.0550
0.3400, 1.0000, 0.3000, 0.0830, 0.1850, 0.1072, 0.5300, 0.3000, 0.4820, 0.0920, 0.0850
0.6700, 0.0000, 0.2070, 0.0830, 0.1700, 0.0998, 0.5900, 0.3000, 0.4025, 0.0770, 0.0890
0.4900, 0.0000, 0.2560, 0.0760, 0.1610, 0.0998, 0.5100, 0.3000, 0.3932, 0.0780, 0.0310
0.5500, 1.0000, 0.2290, 0.0810, 0.1230, 0.0672, 0.4100, 0.3000, 0.4304, 0.0880, 0.1290
0.5900, 1.0000, 0.2510, 0.0900, 0.1630, 0.1014, 0.4600, 0.4000, 0.4357, 0.0910, 0.0830
0.5300, 0.0000, 0.3320, 0.0827, 0.1860, 0.1068, 0.4600, 0.4040, 0.5112, 0.1020, 0.2750
0.4800, 1.0000, 0.2410, 0.1100, 0.2090, 0.1346, 0.5800, 0.4000, 0.4407, 0.1000, 0.0650
0.5200, 0.0000, 0.2950, 0.1043, 0.2110, 0.1328, 0.4900, 0.4310, 0.4984, 0.0980, 0.1980
0.6900, 0.0000, 0.2960, 0.1220, 0.2310, 0.1284, 0.5600, 0.4000, 0.5451, 0.0860, 0.2360
0.6000, 1.0000, 0.2280, 0.1100, 0.2450, 0.1898, 0.3900, 0.6000, 0.4394, 0.0880, 0.2530
0.4600, 1.0000, 0.2270, 0.0830, 0.1830, 0.1258, 0.3200, 0.6000, 0.4836, 0.0750, 0.1240
0.5100, 1.0000, 0.2620, 0.1010, 0.1610, 0.0996, 0.4800, 0.3000, 0.4205, 0.0880, 0.0440
0.6700, 1.0000, 0.2350, 0.0960, 0.2070, 0.1382, 0.4200, 0.5000, 0.4898, 0.1110, 0.1720
0.4900, 0.0000, 0.2210, 0.0850, 0.1360, 0.0634, 0.6200, 0.2190, 0.3970, 0.0720, 0.1140
0.4600, 1.0000, 0.2650, 0.0940, 0.2470, 0.1602, 0.5900, 0.4000, 0.4935, 0.1110, 0.1420
0.4700, 0.0000, 0.3240, 0.1050, 0.1880, 0.1250, 0.4600, 0.4090, 0.4443, 0.0990, 0.1090
0.7500, 0.0000, 0.3010, 0.0780, 0.2220, 0.1542, 0.4400, 0.5050, 0.4779, 0.0970, 0.1800
0.2800, 0.0000, 0.2420, 0.0930, 0.1740, 0.1064, 0.5400, 0.3000, 0.4220, 0.0840, 0.1440
0.6500, 1.0000, 0.3130, 0.1100, 0.2130, 0.1280, 0.4700, 0.5000, 0.5247, 0.0910, 0.1630
0.4200, 0.0000, 0.3010, 0.0910, 0.1820, 0.1148, 0.4900, 0.4000, 0.4511, 0.0820, 0.1470
0.5100, 0.0000, 0.2450, 0.0790, 0.2120, 0.1286, 0.6500, 0.3000, 0.4522, 0.0910, 0.0970
0.5300, 1.0000, 0.2770, 0.0950, 0.1900, 0.1018, 0.4100, 0.5000, 0.5464, 0.1010, 0.2200
0.5400, 0.0000, 0.2320, 0.1107, 0.2380, 0.1628, 0.4800, 0.4960, 0.4913, 0.1080, 0.1900
0.7300, 0.0000, 0.2700, 0.1020, 0.2110, 0.1210, 0.6700, 0.3000, 0.4745, 0.0990, 0.1090
0.5400, 0.0000, 0.2680, 0.1080, 0.1760, 0.0806, 0.6700, 0.3000, 0.4956, 0.1060, 0.1910
0.4200, 0.0000, 0.2920, 0.0930, 0.2490, 0.1742, 0.4500, 0.6000, 0.5004, 0.0920, 0.1220
0.7500, 0.0000, 0.3120, 0.1177, 0.2290, 0.1388, 0.2900, 0.7900, 0.5724, 0.1060, 0.2300
0.5500, 1.0000, 0.3210, 0.1127, 0.2070, 0.0924, 0.2500, 0.8280, 0.6105, 0.1110, 0.2420
0.6800, 1.0000, 0.2570, 0.1090, 0.2330, 0.1126, 0.3500, 0.7000, 0.6057, 0.1050, 0.2480
0.5700, 0.0000, 0.2690, 0.0980, 0.2460, 0.1652, 0.3800, 0.7000, 0.5366, 0.0960, 0.2490
0.4800, 0.0000, 0.3140, 0.0753, 0.2420, 0.1516, 0.3800, 0.6370, 0.5568, 0.1030, 0.1920
0.6100, 1.0000, 0.2560, 0.0850, 0.1840, 0.1162, 0.3900, 0.5000, 0.4970, 0.0980, 0.1310
0.6900, 0.0000, 0.3700, 0.1030, 0.2070, 0.1314, 0.5500, 0.4000, 0.4635, 0.0900, 0.2370
0.3800, 0.0000, 0.3260, 0.0770, 0.1680, 0.1006, 0.4700, 0.4000, 0.4625, 0.0960, 0.0780
0.4500, 1.0000, 0.2120, 0.0940, 0.1690, 0.0968, 0.5500, 0.3000, 0.4454, 0.1020, 0.1350
0.5100, 1.0000, 0.2920, 0.1070, 0.1870, 0.1390, 0.3200, 0.6000, 0.4382, 0.0950, 0.2440
0.7100, 1.0000, 0.2400, 0.0840, 0.1380, 0.0858, 0.3900, 0.4000, 0.4190, 0.0900, 0.1990
0.5700, 0.0000, 0.3610, 0.1170, 0.1810, 0.1082, 0.3400, 0.5000, 0.5268, 0.1000, 0.2700
0.5600, 1.0000, 0.2580, 0.1030, 0.1770, 0.1144, 0.3400, 0.5000, 0.4963, 0.0990, 0.1640
0.3200, 1.0000, 0.2200, 0.0880, 0.1370, 0.0786, 0.4800, 0.3000, 0.3951, 0.0780, 0.0720
0.5000, 0.0000, 0.2190, 0.0910, 0.1900, 0.1112, 0.6700, 0.3000, 0.4078, 0.0770, 0.0960
0.4300, 0.0000, 0.3430, 0.0840, 0.2560, 0.1726, 0.3300, 0.8000, 0.5529, 0.1040, 0.3060
0.5400, 1.0000, 0.2520, 0.1150, 0.1810, 0.1200, 0.3900, 0.5000, 0.4701, 0.0920, 0.0910
0.3100, 0.0000, 0.2330, 0.0850, 0.1900, 0.1308, 0.4300, 0.4000, 0.4394, 0.0770, 0.2140
0.5600, 0.0000, 0.2570, 0.0800, 0.2440, 0.1516, 0.5900, 0.4000, 0.5118, 0.0950, 0.0950
0.4400, 0.0000, 0.2510, 0.1330, 0.1820, 0.1130, 0.5500, 0.3000, 0.4249, 0.0840, 0.2160
0.5700, 1.0000, 0.3190, 0.1110, 0.1730, 0.1162, 0.4100, 0.4000, 0.4369, 0.0870, 0.2630

Test data:


# diabetes_norm_test_100.txt
#
0.6400, 1.0000, 0.2840, 0.1110, 0.1840, 0.1270, 0.4100, 0.4000, 0.4382, 0.0970, 0.1780
0.4300, 0.0000, 0.2810, 0.1210, 0.1920, 0.1210, 0.6000, 0.3000, 0.4007, 0.0930, 0.1130
0.1900, 0.0000, 0.2530, 0.0830, 0.2250, 0.1566, 0.4600, 0.5000, 0.4719, 0.0840, 0.2000
0.7100, 1.0000, 0.2610, 0.0850, 0.2200, 0.1524, 0.4700, 0.5000, 0.4635, 0.0910, 0.1390
0.5000, 1.0000, 0.2800, 0.1040, 0.2820, 0.1968, 0.4400, 0.6000, 0.5328, 0.0950, 0.1390
0.5900, 1.0000, 0.2360, 0.0730, 0.1800, 0.1074, 0.5100, 0.4000, 0.4682, 0.0840, 0.0880
0.5700, 0.0000, 0.2450, 0.0930, 0.1860, 0.0966, 0.7100, 0.3000, 0.4522, 0.0910, 0.1480
0.4900, 1.0000, 0.2100, 0.0820, 0.1190, 0.0854, 0.2300, 0.5000, 0.3970, 0.0740, 0.0880
0.4100, 1.0000, 0.3200, 0.1260, 0.1980, 0.1042, 0.4900, 0.4000, 0.5412, 0.1240, 0.2430
0.2500, 1.0000, 0.2260, 0.0850, 0.1300, 0.0710, 0.4800, 0.3000, 0.4007, 0.0810, 0.0710
0.5200, 1.0000, 0.1970, 0.0810, 0.1520, 0.0534, 0.8200, 0.2000, 0.4419, 0.0820, 0.0770
0.3400, 0.0000, 0.2120, 0.0840, 0.2540, 0.1134, 0.5200, 0.5000, 0.6094, 0.0920, 0.1090
0.4200, 1.0000, 0.3060, 0.1010, 0.2690, 0.1722, 0.5000, 0.5000, 0.5455, 0.1060, 0.2720
0.2800, 1.0000, 0.2550, 0.0990, 0.1620, 0.1016, 0.4600, 0.4000, 0.4277, 0.0940, 0.0600
0.4700, 1.0000, 0.2330, 0.0900, 0.1950, 0.1258, 0.5400, 0.4000, 0.4331, 0.0730, 0.0540
0.3200, 1.0000, 0.3100, 0.1000, 0.1770, 0.0962, 0.4500, 0.4000, 0.5187, 0.0770, 0.2210
0.4300, 0.0000, 0.1850, 0.0870, 0.1630, 0.0936, 0.6100, 0.2670, 0.3738, 0.0800, 0.0900
0.5900, 1.0000, 0.2690, 0.1040, 0.1940, 0.1266, 0.4300, 0.5000, 0.4804, 0.1060, 0.3110
0.5300, 0.0000, 0.2830, 0.1010, 0.1790, 0.1070, 0.4800, 0.4000, 0.4788, 0.1010, 0.2810
0.6000, 0.0000, 0.2570, 0.1030, 0.1580, 0.0846, 0.6400, 0.2000, 0.3850, 0.0970, 0.1820
0.5400, 1.0000, 0.3610, 0.1150, 0.1630, 0.0984, 0.4300, 0.4000, 0.4682, 0.1010, 0.3210
0.3500, 1.0000, 0.2410, 0.0947, 0.1550, 0.0974, 0.3200, 0.4840, 0.4852, 0.0940, 0.0580
0.4900, 1.0000, 0.2580, 0.0890, 0.1820, 0.1186, 0.3900, 0.5000, 0.4804, 0.1150, 0.2620
0.5800, 0.0000, 0.2280, 0.0910, 0.1960, 0.1188, 0.4800, 0.4000, 0.4984, 0.1150, 0.2060
0.3600, 1.0000, 0.3910, 0.0900, 0.2190, 0.1358, 0.3800, 0.6000, 0.5421, 0.1030, 0.2330
0.4600, 1.0000, 0.4220, 0.0990, 0.2110, 0.1370, 0.4400, 0.5000, 0.5011, 0.0990, 0.2420
0.4400, 1.0000, 0.2660, 0.0990, 0.2050, 0.1090, 0.4300, 0.5000, 0.5580, 0.1110, 0.1230
0.4600, 0.0000, 0.2990, 0.0830, 0.1710, 0.1130, 0.3800, 0.4500, 0.4585, 0.0980, 0.1670
0.5400, 0.0000, 0.2100, 0.0780, 0.1880, 0.1074, 0.7000, 0.3000, 0.3970, 0.0730, 0.0630
0.6300, 1.0000, 0.2550, 0.1090, 0.2260, 0.1032, 0.4600, 0.5000, 0.5951, 0.0870, 0.1970
0.4100, 1.0000, 0.2420, 0.0900, 0.1990, 0.1236, 0.5700, 0.4000, 0.4522, 0.0860, 0.0710
0.2800, 0.0000, 0.2540, 0.0930, 0.1410, 0.0790, 0.4900, 0.3000, 0.4174, 0.0910, 0.1680
0.1900, 0.0000, 0.2320, 0.0750, 0.1430, 0.0704, 0.5200, 0.3000, 0.4635, 0.0720, 0.1400
0.6100, 1.0000, 0.2610, 0.1260, 0.2150, 0.1298, 0.5700, 0.4000, 0.4949, 0.0960, 0.2170
0.4800, 0.0000, 0.3270, 0.0930, 0.2760, 0.1986, 0.4300, 0.6420, 0.5148, 0.0910, 0.1210
0.5400, 1.0000, 0.2730, 0.1000, 0.2000, 0.1440, 0.3300, 0.6000, 0.4745, 0.0760, 0.2350
0.5300, 1.0000, 0.2660, 0.0930, 0.1850, 0.1224, 0.3600, 0.5000, 0.4890, 0.0820, 0.2450
0.4800, 0.0000, 0.2280, 0.1010, 0.1100, 0.0416, 0.5600, 0.2000, 0.4127, 0.0970, 0.0400
0.5300, 0.0000, 0.2880, 0.1117, 0.1450, 0.0872, 0.4600, 0.3150, 0.4078, 0.0850, 0.0520
0.2900, 1.0000, 0.1810, 0.0730, 0.1580, 0.0990, 0.4100, 0.4000, 0.4500, 0.0780, 0.1040
0.6200, 0.0000, 0.3200, 0.0880, 0.1720, 0.0690, 0.3800, 0.4000, 0.5784, 0.1000, 0.1320
0.5000, 1.0000, 0.2370, 0.0920, 0.1660, 0.0970, 0.5200, 0.3000, 0.4443, 0.0930, 0.0880
0.5800, 1.0000, 0.2360, 0.0960, 0.2570, 0.1710, 0.5900, 0.4000, 0.4905, 0.0820, 0.0690
0.5500, 1.0000, 0.2460, 0.1090, 0.1430, 0.0764, 0.5100, 0.3000, 0.4357, 0.0880, 0.2190
0.5400, 0.0000, 0.2260, 0.0900, 0.1830, 0.1042, 0.6400, 0.3000, 0.4304, 0.0920, 0.0720
0.3600, 0.0000, 0.2780, 0.0730, 0.1530, 0.1044, 0.4200, 0.4000, 0.3497, 0.0730, 0.2010
0.6300, 1.0000, 0.2410, 0.1110, 0.1840, 0.1122, 0.4400, 0.4000, 0.4935, 0.0820, 0.1100
0.4700, 1.0000, 0.2650, 0.0700, 0.1810, 0.1048, 0.6300, 0.3000, 0.4190, 0.0700, 0.0510
0.5100, 1.0000, 0.3280, 0.1120, 0.2020, 0.1006, 0.3700, 0.5000, 0.5775, 0.1090, 0.2770
0.4200, 0.0000, 0.1990, 0.0760, 0.1460, 0.0832, 0.5500, 0.3000, 0.3664, 0.0790, 0.0630
0.3700, 1.0000, 0.2360, 0.0940, 0.2050, 0.1388, 0.5300, 0.4000, 0.4190, 0.1070, 0.1180
0.2800, 0.0000, 0.2210, 0.0820, 0.1680, 0.1006, 0.5400, 0.3000, 0.4205, 0.0860, 0.0690
0.5800, 0.0000, 0.2810, 0.1110, 0.1980, 0.0806, 0.3100, 0.6000, 0.6068, 0.0930, 0.2730
0.3200, 0.0000, 0.2650, 0.0860, 0.1840, 0.1016, 0.5300, 0.4000, 0.4990, 0.0780, 0.2580
0.2500, 1.0000, 0.2350, 0.0880, 0.1430, 0.0808, 0.5500, 0.3000, 0.3584, 0.0830, 0.0430
0.6300, 0.0000, 0.2600, 0.0857, 0.1550, 0.0782, 0.4600, 0.3370, 0.5037, 0.0970, 0.1980
0.5200, 0.0000, 0.2780, 0.0850, 0.2190, 0.1360, 0.4900, 0.4000, 0.5136, 0.0750, 0.2420
0.6500, 1.0000, 0.2850, 0.1090, 0.2010, 0.1230, 0.4600, 0.4000, 0.5075, 0.0960, 0.2320
0.4200, 0.0000, 0.3060, 0.1210, 0.1760, 0.0928, 0.6900, 0.3000, 0.4263, 0.0890, 0.1750
0.5300, 0.0000, 0.2220, 0.0780, 0.1640, 0.0810, 0.7000, 0.2000, 0.4174, 0.1010, 0.0930
0.7900, 1.0000, 0.2330, 0.0880, 0.1860, 0.1284, 0.3300, 0.6000, 0.4812, 0.1020, 0.1680
0.4300, 0.0000, 0.3540, 0.0930, 0.1850, 0.1002, 0.4400, 0.4000, 0.5318, 0.1010, 0.2750
0.4400, 0.0000, 0.3140, 0.1150, 0.1650, 0.0976, 0.5200, 0.3000, 0.4344, 0.0890, 0.2930
0.6200, 1.0000, 0.3780, 0.1190, 0.1130, 0.0510, 0.3100, 0.4000, 0.5043, 0.0840, 0.2810
0.3300, 0.0000, 0.1890, 0.0700, 0.1620, 0.0918, 0.5900, 0.3000, 0.4025, 0.0580, 0.0720
0.5600, 0.0000, 0.3500, 0.0793, 0.1950, 0.1408, 0.4200, 0.4640, 0.4111, 0.0960, 0.1400
0.6600, 0.0000, 0.2170, 0.1260, 0.2120, 0.1278, 0.4500, 0.4710, 0.5278, 0.1010, 0.1890
0.3400, 1.0000, 0.2530, 0.1110, 0.2300, 0.1620, 0.3900, 0.6000, 0.4977, 0.0900, 0.1810
0.4600, 1.0000, 0.2380, 0.0970, 0.2240, 0.1392, 0.4200, 0.5000, 0.5366, 0.0810, 0.2090
0.5000, 0.0000, 0.3180, 0.0820, 0.1360, 0.0692, 0.5500, 0.2000, 0.4078, 0.0850, 0.1360
0.6900, 0.0000, 0.3430, 0.1130, 0.2000, 0.1238, 0.5400, 0.4000, 0.4710, 0.1120, 0.2610
0.3400, 0.0000, 0.2630, 0.0870, 0.1970, 0.1200, 0.6300, 0.3000, 0.4249, 0.0960, 0.1130
0.7100, 1.0000, 0.2700, 0.0933, 0.2690, 0.1902, 0.4100, 0.6560, 0.5242, 0.0930, 0.1310
0.4700, 0.0000, 0.2720, 0.0800, 0.2080, 0.1456, 0.3800, 0.6000, 0.4804, 0.0920, 0.1740
0.4100, 0.0000, 0.3380, 0.1233, 0.1870, 0.1270, 0.4500, 0.4160, 0.4318, 0.1000, 0.2570
0.3400, 0.0000, 0.3300, 0.0730, 0.1780, 0.1146, 0.5100, 0.3490, 0.4127, 0.0920, 0.0550
0.5100, 0.0000, 0.2410, 0.0870, 0.2610, 0.1756, 0.6900, 0.4000, 0.4407, 0.0930, 0.0840
0.4300, 0.0000, 0.2130, 0.0790, 0.1410, 0.0788, 0.5300, 0.3000, 0.3829, 0.0900, 0.0420
0.5500, 0.0000, 0.2300, 0.0947, 0.1900, 0.1376, 0.3800, 0.5000, 0.4277, 0.1060, 0.1460
0.5900, 1.0000, 0.2790, 0.1010, 0.2180, 0.1442, 0.3800, 0.6000, 0.5187, 0.0950, 0.2120
0.2700, 1.0000, 0.3360, 0.1100, 0.2460, 0.1566, 0.5700, 0.4000, 0.5088, 0.0890, 0.2330
0.5100, 1.0000, 0.2270, 0.1030, 0.2170, 0.1624, 0.3000, 0.7000, 0.4812, 0.0800, 0.0910
0.4900, 1.0000, 0.2740, 0.0890, 0.1770, 0.1130, 0.3700, 0.5000, 0.4905, 0.0970, 0.1110
0.2700, 0.0000, 0.2260, 0.0710, 0.1160, 0.0434, 0.5600, 0.2000, 0.4419, 0.0790, 0.1520
0.5700, 1.0000, 0.2320, 0.1073, 0.2310, 0.1594, 0.4100, 0.5630, 0.5030, 0.1120, 0.1200
0.3900, 1.0000, 0.2690, 0.0930, 0.1360, 0.0754, 0.4800, 0.3000, 0.4143, 0.0990, 0.0670
0.6200, 1.0000, 0.3460, 0.1200, 0.2150, 0.1292, 0.4300, 0.5000, 0.5366, 0.1230, 0.3100
0.3700, 0.0000, 0.2330, 0.0880, 0.2230, 0.1420, 0.6500, 0.3400, 0.4357, 0.0820, 0.0940
0.4600, 0.0000, 0.2110, 0.0800, 0.2050, 0.1444, 0.4200, 0.5000, 0.4533, 0.0870, 0.1830
0.6800, 1.0000, 0.2350, 0.1010, 0.1620, 0.0854, 0.5900, 0.3000, 0.4477, 0.0910, 0.0660
0.5100, 0.0000, 0.3150, 0.0930, 0.2310, 0.1440, 0.4900, 0.4700, 0.5252, 0.1170, 0.1730
0.4100, 0.0000, 0.2080, 0.0860, 0.2230, 0.1282, 0.8300, 0.3000, 0.4078, 0.0890, 0.0720
0.5300, 0.0000, 0.2650, 0.0970, 0.1930, 0.1224, 0.5800, 0.3000, 0.4143, 0.0990, 0.0490
0.4500, 0.0000, 0.2420, 0.0830, 0.1770, 0.1184, 0.4500, 0.4000, 0.4220, 0.0820, 0.0640
0.3300, 0.0000, 0.1950, 0.0800, 0.1710, 0.0854, 0.7500, 0.2000, 0.3970, 0.0800, 0.0480
0.6000, 1.0000, 0.2820, 0.1120, 0.1850, 0.1138, 0.4200, 0.4000, 0.4984, 0.0930, 0.1780
0.4700, 1.0000, 0.2490, 0.0750, 0.2250, 0.1660, 0.4200, 0.5000, 0.4443, 0.1020, 0.1040
0.6000, 1.0000, 0.2490, 0.0997, 0.1620, 0.1066, 0.4300, 0.3770, 0.4127, 0.0950, 0.1320
0.3600, 0.0000, 0.3000, 0.0950, 0.2010, 0.1252, 0.4200, 0.4790, 0.5130, 0.0850, 0.2200
0.3600, 0.0000, 0.1960, 0.0710, 0.2500, 0.1332, 0.9700, 0.3000, 0.4595, 0.0920, 0.0570
Posted in Machine Learning | Leave a comment

Ridge Regression From Scratch Using Python With SGD Training

One morning before work, I noticed that it had been several months since I last implemented linear ridge regression, from scratch, using Python. I figured I’d do so, using the style of the scikit-learn Ridge module.

Ridge regression (sometimes called linear ridge regression, but is much different from kernel ridge regression) is just basic linear regression with L2 regularization. L2 regularization is applied during training and discourages model weights from becoming extremely large, because large model weights indicate likely model overfitting.

Mathematically, L2 regularization adds a penalty equal to the sum of the squared weights in the model. When you compute a gradient for stochastic gradient descent training, the squared part goes away and you end up adding alpha * weight[j] to the gradient for each weight. The alpha is a small constant, typically about 0.01, that moderates the amount of regularization.


Details

If you use closed form training — left pseudo-inverse via normal equations with Cholesky, or Moore-Penrose pseudo-inverse with SVD — then L2 is accomplished by adding alpha to the diagonal elements of the primary matrix before applying the matrix inverse operation — a completely different approach.

In math literature, the L2 constant is called lambda, but in computer science we use alpha to avoid colliding with the lambda keyword that is found in many programming languages.

There is also L1 regularization, but I never use it — explaining why is a non-trial topic.

In math notation, the number of rows in a matrix is usually called m and the number of columns is n. I use non-standard notation of n for the number of rows, and dim for the number of columns.

My demo implementation uses explicit for-loops for clarity. Numpy supports all kinds of bulk operations on vectors and matrices, which is more efficient, but much less clear.

Years ago, the error term was sometimes computed as pred_y – target_y and sometimes computed as target_y – pred_y. The pred_y – target_y form seems universal now. If you use target_y – pred_y, then the -= in the weight update becomes += (which always seemed a bit more logical to me, but I sadly was never elected King of All Machine Learning).

The scikit Ridge module uses a clever mini-algorithm to automatically estimate a good value for the learning rate, but I suspect that the clever mini-algorithm algorithm works well only with the weird SAG training algorithm. You can always use a simple grid search to find a good learning rate.

Both the scikit Ridge module and my from-scratch version use a very simple early-exit for the SGD training loop. The check is if the smallest change in one of the weights, divided by the smallest one of the current weights, is less than some tolerance, then early exit. This technique is simple but risks a too-early exit. A more robust approach is to only early exit if the small-change criterion is met 3 (or perhaps 5) consecutive times.

I use A[i][j] matrix syntax rather than the slightly more efficient A[i,j] syntax — just a personal preference.

To evaluate the from-scratch ridge regression model and the scikit Ridge model, I implemented a program-defined mse() function and an accuracy() function. I allowed for the fact that all scikit model predict() methods require a matrix as input and return a vector as output, rather than accepting a vector as input and returning a single value as output.


The key training code looks like:

for-each epoch (pass through all data)
  for-each training item (in scrambled order)
    get input x
    get target y
    compute pred y using curr weights
    err = pred y - target y
    for-each weight[j]
      l2_term = alpha * weight[j]
      gradient = err * x[j] + l2_term
      weights[j] -= lrn_rate * gradient
    end-loop
    bias -= lrn_rate * err (no L2 for bias)
  end-loop
end-loop

For my demo, I used a set of synthetic data that 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 first five values on each line are predictors. The last value on each line is the target y value to predict. The data has been pre-scaled so that all predictors are between -1 and +1. When using ridge regression, you should scale data so that all weights are penalized equally. The scikit StandardScaler module can do this if your data is not pre-scaled.

There are 200 training items and 40 test items. The data has complex non-linear relations so linear regression, ridge or plain, will not work very well.

The key parts of the demo output are:

Begin SGD ridge regression scratch Python

Loading train (200) and test (40) data

Creating scratch Python ridge regression model
Done

Training model with SGD
Setting lrn_rate = 0.0010
Setting max_epochs = 100
Setting L2 alpha = 0.005000
Setting early-exit tol = 0.000100

epoch:     0   MSE =   0.1133
epoch:    20   MSE =   0.0045
epoch:    40   MSE =   0.0027
epoch:    60   MSE =   0.0026
epoch:    80   MSE =   0.0026
Done. Used 94 epochs

Model weights:
[-0.2614  0.0331 -0.0457  0.0351 -0.1129]
Model bias = 0.3618

MSE train = 0.00259
MSE test = 0.00193

The learning rate, max_epochs, L2 alpha, and early training exit tolerance parameters all had to be determined by trial and error.

To sort-of validate my demo system, I ran the data through the scikit Ridge module. The model produced by scikit was almost identical to the from-scratch model:

Model weights:
[-0.2618  0.0331 -0.0453  0.0353 -0.1132]
Model bias = 0.3618

MSE train = 0.00259
MSE test = 0.00194

One weird thing about the scikit Ridge module is that it doesn’t use alpha directly. Instead it scales alpha by dividing it by the number of training items. So, when using scikit, if you set alpha = 1.0 (the default) with 200 training items, the alpha value used behind the scenes is 1.0 / 200 = 0.005. This is somewhat confusing, because the scikit documentation doesn’t explain this (I had to wade through the scikit source code to figure this out).

The scikit Ridge module can use one of seven training algorithms: ‘auto’, ‘svd’, ‘cholesky’, ‘lsqr’, ‘sparse_cg’, ‘sag’, ‘saga, ‘lbfgs’. The ‘auto’ means use some mini-algorithm to guess the best training algorithm. My demo specified ‘sag’ (stochastic average gradient) because it is a strange variation of SGD that is optimized for speed.



Machine learning linear regresssion techniques were developed in the late 1950s and early 1960s. The space race to the moon was also going on in those years. It was a thrilling time to be alive — I was there.

Left: This illustration is from the June 27, 1953 issue of Collier’s Magazine, as part of a multi-issue series “Man Will Conquer Space Soon!” The design is loosely based on one by famous rocket scientist Werner von Braun. The artist is Chesley Bonstell. I give the illustration my personal A- grade.

Center and Right: The movie “War of the Satellites” (1958) almost certainly used the Collier’s illustration as the inspiration for its special effects. Aliens who are just energy (no bodies) plan to invade Earth. They don’t succeed. A very low-budget film but one that has a certain charm (for me anyway). Most critics give this movie the equivalent of a D grade, but I inexplicably give the movie my personal C+ grade. I am easily entertained.


Demo program. Replace the “lt” in the accuracy() function with the Boolean less-than operator.

# ridge_regression_sgd_scratch.py
# Stochastic Gradient Descent training

import numpy as np

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

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

    if scikit == False:
      pred_y = model.predict(x)
    else:
      pred_y = model.predict([x])[0]

    if np.abs(pred_y - 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, scikit=False):
  n = len(data_X)
  sum = 0.0
  for i in range(n):
    x = data_X[i]
    y = data_y[i]
    if scikit == False:
      pred_y = model.predict(x)
    else:
      pred_y = model.predict([x])[0]
    diff = pred_y - y
    sum += diff * diff
  return sum /n

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

class RidgeRegressor:
  # ridge regression using SGD
  def __init__(self, seed=0):
    self.weights = None
    self.bias = 0.0
    self.rnd = np.random.RandomState(seed)

  def predict(self, x):
    # x is a vector
    sum = 0.0
    for j in range(len(x)):
      sum += self.weights[j] * x[j]
    sum += self.bias
    return sum

  def train(self, train_X, train_y, lrn_rate=0.001,
    max_epochs=1000, tol=0.0001, alpha=0.01):
    # non-scikit style: alpha is not scaled
    n = len(train_X)
    dim = len(train_X[0])

    self.weights = np.zeros(dim) # init wts to small rnd
    lo = -0.01; hi = 0.01
    for j in range(dim):
      self.weights[j] = (hi - lo) * self.rnd.random() + lo

    prev_weights = np.zeros(dim)
    indices = np.arange(n)

    for epoch in range(max_epochs): # each epoch
      self.rnd.shuffle(indices)
      prev_weights = np.copy(self.weights)
      for i in range(n): # each item
        idx = indices[i]
        x = train_X[idx]
        y = train_y[idx]
        pred_y = self.predict(x)
        err = pred_y - y  # most common form
        
        for j in range(dim): # each weight
          l2_term = alpha * self.weights[j]
          self.weights[j] -= lrn_rate * (err * x[j] + l2_term)
        self.bias -= lrn_rate * err  # no L2 for bias

      # display progress 5 times
      if epoch % (max_epochs // 5) == 0:
        mse = self.mse(train_X, train_y)
        s1 = "epoch: %5d" % epoch
        s2 = "   MSE = %8.4f" % mse
        print(s1 + s2) 

      # check for early stop
      max_change_in_wts = np.abs(self.weights - \
        prev_weights).max()
      max_weight = np.abs(self.weights).max()
      if max_change_in_wts / max_weight < tol:
        return epoch
      
    return epoch

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

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

print("\nBegin SGD ridge regression scratch Python ")

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

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

print("\n========================== ")

print("\nCreating scratch Python ridge regression model ")
model = RidgeRegressor()
print("Done ")

print("\nTraining model with SGD ")
lrn_rate = 0.001
max_epochs = 100
alpha = 0.005 # 1/n_items
tol = 0.0001 # for early-stop
print("Setting lrn_rate = %0.4f " % lrn_rate)
print("Setting max_epochs = " + str(max_epochs))
print("Setting L2 alpha = %0.6f " % alpha)
print("Setting early-exit tol = %0.6f \n" % tol)
epochs = model.train(train_X, train_y, lrn_rate, 
  max_epochs, tol, alpha)
print("Done. Used " + str(epochs) + " epochs " )

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

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.5f " % mse_train)
print("MSE test = %0.5f " % mse_test)

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

print("\n========================== ")

from sklearn.linear_model import Ridge

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

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

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

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

acc_train = \
  accuracy(model, train_X, train_y, 0.10, scikit=True)
acc_test = \
  accuracy(model, test_X, test_y, 0.10, scikit=True)
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, scikit=True)
mse_test = mse(model, test_X, test_y, scikit=True)
print("\nMSE train = %0.5f " % mse_train)
print("MSE test = %0.5f " % mse_test)

print("\nEnd ridge regression 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 | 1 Comment

Matrix Left Pseudo-Inverse Via Normal Equations (Cholesky) – Refactor and Test Using C#

In machine learning, to train a linear model (linear regression, quadratic regression, etc.) you can use several techniques including stochastic gradient descent (SGD), L-BFGS optimization, and pseudo-inverse. There are two main categories of pseudo-inverse algorithms for machine learning scenarios where the training data has more rows than columns: relaxed Moore-Penrose pseudo-inverse, and left pseudo-inverse via normal equations.

For relaxed Moore-Penrose pseudo-inverse there are two main categories of techniques: using singular value decomposition (SVD) or QR decomposition. There are several SVD algorithms, including Jacobi, Householder+QR, Golub-Kahan-Reinsch, and Lanczos. There are several QR algorithms including Householder, modified Gram-Schmidt, and Givens.

For left pseudo-inverse via normal equations, the most common technique is Cholesky decomposition. There are several versions of Cholesky, but the most common is the Banachiewicz algorithm.

Each of these pseudo-inverse techniques has pros and cons related to numerical stability for large matrices, performance speed, and implementation complexity. For example, SVD Jacobi is the most stable (in theory) but the most complex. And QR modified Gram-Schmidt is relatively simple but less stable (in theory). And Cholesky is very stable but is vulnerable to poorly conditioned matrices.

I write code almost every day. Coding is a skill that must be practiced, and I enjoy coding. One morning before work, I figured I’d refactor my C# implementation of left pseudo-inverse via normal equations using Cholesky (Banachiewicz). So I did.

If you have a matrix A of training data (or a design matrix derived from the training data) with more rows than columns, the left pseudo-inverse can be computed as pinv = inv(At * A) * At, where At is the transpose of A, * is matrix multiplication, and inv() is any matrix inverse function. But At * A is square, symmetric, positive definite (assuming you can a small conditioning constant to the diagonal elements) and you can use a relatively simple matrix inverse — Cholesky decomposition.

This technique — left pseudo-inverse via normal equations using Cholesky — is dramatically simpler than any other form of pseudo-inverse for machine learning scenarios. But the technique can fail when the source matrix has many very small values or many large values. The problem is the At * A computation, which has lots of multiplication operations.

The key code in my demo is:

public class Cholesky
{
  public static double[][] MatPseudoInv(double[][] A)
  {
    // pinv = inv(At * A) * At
    double[][] At = MatTranspose(A);
    double[][] AtA = MatProduct(At, A);
    double[][] AtAinv = MatInvCholesky(AtA);
    double[][] result = MatProduct(AtAinv, At);
    return result;
  }

  private static double[][] MatInvCholesky(double[][] A)
  {
    // A must be square symmetric positive definite
    // A = L * Lt, inv(A) = inv(Lt) * inv(L)
    double[][] L = MatDecompCholesky(A); // lower tri
    double[][] Lt = MatTranspose(L); // upper tri
    double[][] invL = MatInvLowerTri(L);
    double[][] invLt = MatInvUpperTri(Lt);
    double[][] result = MatProduct(invLt, invL);
    return result;
  }

  private static double[][] MatDecompCholesky(double[][] A)
  {
    // lots of code
  }

  // lots of minor helpers like MatTranspose()

} // class Cholesky

I tested my refactored pseudo-inverse implementation by generating 10,000 random matrices, computing the pseudo-inverse, then verifying that A * Apinv * A = A to within 1.0e-8 in each cell. Each matrix has between 100 and 1,000 rows, between 2 and 20 columns, and each value is between -10.0 and +10.0. All 10,000 tests passed. This is by no means thorough testing, but it does indicate that the implementation is basically correct for non-pathological matrices.

Good fun.



I miss the days of newspapers.


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

using System.IO;

namespace MatrixPseudoInverseCholesky
{
  internal class MatrixPseudoInverseCholeskyProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin left pseudo-inverse" +
        " via normal equations (Cholesky) ");

      int minRows = 100; int maxRows = 1000;
      int minCols = 2; int maxCols = 20;
      Random rnd = new Random(0);
      int nTrials = 10000;
      int nPass = 0; int nFail = 0;
      for (int i = 0; i "lt" nTrials; ++i)
      {
        Console.WriteLine("\n==========");
        Console.WriteLine("trial " + i);
        double[][] A =
          MakeRandomMatrix(minRows, maxRows,
          minCols, maxCols, rnd);
        double[][] Apinv = Cholesky.MatPseudoInv(A);

        // check (A * Apinv) * A = A
        double[][] AApinv = MatProduct(A, Apinv);
        double[][] C = MatProduct(AApinv, A);

        //Console.WriteLine("\nA = ");
        //MatShow(A, 4, 9);
        //Console.WriteLine("\nApinv = ");
        //MatShow(Apinv, 4, 9);
        //Console.WriteLine("\n(A * Apinv) * A = ");
        //MatShow(C, 4, 9);

        if (MatAreEqual(C, A, 1.0e-8) == true)
        {
          Console.WriteLine("pass");
          ++nPass;
        }
        else
        {
          Console.WriteLine("FAIL");
          Console.WriteLine("nRows = " + A.Length +
            " nCols = " + A[0].Length);
          Console.ReadLine();
          ++nFail;
        }
        Console.WriteLine("==========");
     } // nTrials

      Console.WriteLine("\nNumber pass = " + nPass);
      Console.WriteLine("Number fail = " + nFail);

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

    // ------------------------------------------------------
    // helpers: MakeRandomMatrix(), MatProduct(), 
    // MatShow(), MatAreEqual()
    // ------------------------------------------------------

    public static double[][] MakeRandomMatrix(int minRows,
      int maxRows, int minCols, int maxCols, Random rnd)
    {
      int nRows = rnd.Next(minRows, maxRows);  // [2,maxN)
      int nCols = rnd.Next(minCols, maxCols);
      double[][] result = new double[nRows][];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = new double[nCols];

      double lo = -10.0; double hi = 10.0;
      for (int r = 0; r "lt" nRows; ++r)
        for (int c = 0; c "lt" nCols; ++c)
          result[r][c] = (hi - lo) * rnd.NextDouble() + lo;
      return result;
    }

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

    public static double[][] MatProduct(double[][] A,
      double[][] B)
    {
      int aRows = A.Length; int aCols = A[0].Length;
      int bRows = B.Length; int bCols = B[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

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

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

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

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

    public static bool MatAreEqual(double[][] A, double[][] B,
      double epsilon)
    {
      int nRows = A.Length; int nCols = A[0].Length;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          if (Math.Abs(A[i][j] - B[i][j]) "gt" epsilon)
            return false;
      return true;
    }

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

  } // Program

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

  public class Cholesky
  {
    public static double[][] MatPseudoInv(double[][] A)
    {
      // pinv = inv(At * A) * At
      double[][] At = MatTranspose(A);
      double[][] AtA = MatProduct(At, A);
      double[][] AtAinv = MatInvCholesky(AtA);
      double[][] result = MatProduct(AtAinv, At);
      return result;
    }

    private static double[][] MatInvCholesky(double[][] A)
    {
      // A must be square symmetric positive definite
      // A = L * Lt, inv(A) = inv(Lt) * inv(L)
      double[][] L = MatDecompCholesky(A); // lower tri
      double[][] Lt = MatTranspose(L); // upper tri
      double[][] invL = MatInvLowerTri(L);
      double[][] invLt = MatInvUpperTri(Lt);
      double[][] result = MatProduct(invLt, invL);
      return result;
    }

    private static double[][] MatDecompCholesky(double[][] A)
    {
      // A must be square, symmetric, positive definite
      int m = A.Length; int n = A[0].Length;  // m == n
      double[][] L = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        L[i] = new double[n];

      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lte" i; ++j)
        {
          double sum = 0.0;
          for (int k = 0; k "lt" j; ++k)
            sum += L[i][k] * L[j][k];
          if (i == j)
          {
            double tmp = A[i][i] - sum;
            if (tmp "lt" 0.0)
              throw new
                Exception("decomp Cholesky fatal");
            L[i][j] = Math.Sqrt(tmp);
          }
          else
          {
            if (L[j][j] == 0.0)
              throw new
                Exception("decomp Cholesky fatal ");
            L[i][j] = (A[i][j] - sum) / L[j][j];
          }
        } // j
      } // i

      return L;
    }

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

    private static double[][] MatInvLowerTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[k][j] -= result[i][j] * A[k][i];
          }
          result[k][j] /= (A[k][k] + 1.0e-8);
        }
      }
      return result;
    }

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

    private static double[][] MatInvUpperTri(double[][] A)
    {
      int n = A.Length;  // must be square matrix
      double[][] result = MatIdentity(n);

      for (int k = 0; k "lt" n; ++k)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          for (int i = 0; i "lt" k; ++i)
          {
            result[j][k] -= result[j][i] * A[i][k];
          }
          result[j][k] /= (A[k][k] + 1.0e-8); // avoid 0
        }
      }
      return result;
    }

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

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

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

    private static double[][] MatIdentity(int n)
    {
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = 1.0;
      return result;
    }

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

    private static double[][] MatTranspose(double[][] M)
    {
      int nRows = M.Length; int nCols = M[0].Length;
      double[][] result = MatMake(nCols, nRows);  // note
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[j][i] = M[i][j];
      return result;
    }

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

    private static double[][] MatProduct(double[][] A,
      double[][] B)
    {
      int aRows = A.Length; int aCols = A[0].Length;
      int bRows = B.Length; int bCols = B[0].Length;
      if (aCols != bRows)
        throw new Exception("Non-conformable matrices");

      double[][] result = MatMake(aRows, bCols);

      for (int i = 0; i "lt" aRows; ++i) // each row of A
        for (int j = 0; j "lt" bCols; ++j) // each col of B
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

  } // class Cholesky

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

} // ns
Posted in Machine Learning | Leave a comment

I Play Four Interesting Chess Games in a Tournament

When I was in high school (Servite High School in Anaheim, California), I was on the school chess team. Our team — Bob Smith, Tom Law, Tom Quackenbush, Dan Musser, Ed Hernandez, Dan Musser, Dennis Michel, and Mike Ventriglia — was very strong and we won the Orange County championship in our junior and senior years — arguably the de facto State Championship.

We had four players with near-Expert rating strength and a couple with probably-A rating strength. Other schools sometimes had one Master strength player, but those schools couldn’t beat us in a five-board match. For example, Marina High School (Huntington Beach) had Kim Commons (1951-2015), who was a very strong Master.

But then came college, and grad school, and work, and family, and more work — and there was no time for chess for the next several decades.

Well, now I have time and so I played in a local chess tournament on a recent weekend. Four of my games had an interesting tactical sequence, and one game in particular had the longest forced checkmate I’ve ever played — forced mate in 11 moves.

I have mixed feelings about chess. For young people, it teaches planning, risk assessment, mental calculation, sportsmanship, and other good things. But it requires many hundreds of hours of practice to advance beyond beginner level, and some of that time might be spent better on other tasks. And for some people, chess can become a destructive obsession. But the same can be said for almost any activity — video gaming, social media browsing, physical exercise, and on and on.


Black to Move and Mate in 11. In this game, I played the black pieces. My opponent played the London System opening. He attacked on the king side, while I attacked on the queen side. I won with 27… Rxb2! threatening Rc2+. If 28. Kxb2 Rb8+ followed by Qa3 forces checkmate.


Black to Move and Mate in 9. I played the black side of a Sicilian Defense. My opponent played the Snyder Variation (2. b3). My opponent castled queen side, but my queen side attack was quicker than his king side attack. I won with 23… Nxe2+ unleashing the bishop on f6. After 24. Rdxe2 Qxb2+ 25. Kd2 Rd8+ the white king was flushed out to the center of the board and was mated a few moves later.


White to Move and Win the Black Knight. I played the white pieces and opened with 1. e4. My opponent played the rock-solid Caro-Kann Defense. I opted for the Exchange Variation. I won with 35. Ne4 attacking the black queen and eyeing the f6 hole. After 35… Qb5 36. Rxd5 Rxd5 37. Nf6+ Kf8 38. Nxd5 I had won a piece. Knights are tricky in the endgame.


White to Move and Win the Exchange. I played another white side of a Caro-Kann Defense, Exchange Variation. Black was probing weaknesses in my queen side with his knight, but after 19. b4 Nc4 20. Nc6 Qc2 21. Nxb1 I had won a rook for a knight. I eventually stabilized the queen side and broke through to the black king with f4 and f5 but it was a hard fight.


Posted in Miscellaneous | Leave a comment

Two Reasons Why Drop-First Encoding for Linear Regression Categorical Predictor Variables is Preferable to One-Hot Encoding

Bottom line: When using linear regression on data that has categorical predictor variables, you must use drop-first encoding instead of one-hot encoding if you train the model using a closed form pseudo-inverse technique (left pseudo-inverse, Moore-Penrose pseudo-inverse). If you use iterative techniques SGD or L-BFGS to train, you can use one-hot encoding or drop-first encoding, but drop-first is preferred because it gives a more interpretable model than one-hot encoding. In short, for linear regression, drop-first encoding of categorical predictors is the preferred encoding technique.

The goal of a machine learning regression model is to predict a single numeric value. For example, suppose you have data that looks like:

F, 24, michigan, 29500.00, liberal
M, 39, oklahoma, 51200.00, moderate
F, 63, nebraska, 75800.00, conservative
M, 36, virginia, 44500.00, moderate
F, 27, nebraska, 28600.00, liberal
. . .

The fields are sex, age, State (Michigan or Nebraska or Oklahoma or Virginia), annual income, and political leaning. And suppose you want to predict income from the other four variables.

The most fundamental regression technique is linear regression.

In most scenarios, the best way to normalize and encode the data for linear regression is to use divide-by-constant normalization and drop-first encoding:

1, 0.24, 0 0, 0 0 0, 0.29500, 0 1
0, 0.39, 0 1, 0 1 0, 0.51200, 1 0
1, 0.63, 1 0, 1 0 0, 0.75800, 0 0
0, 0.36, 0 1, 0 0 1, 0.44500, 1 0
1, 0.27, 0 0, 1 0 0, 0.28600, 0 1
. . .

The age values are divided by 100, and the income values are divided by 100,000. The binary sex values are encoded as M = 0, F = 1.

If the four State values were one-hot encoded, they would be Michigan = 1000, Nebraska = 0100, Oklahoma = 0010, Virginia = 0001, but drop-first gives Michigan = 000, Nebraska = 100, Oklahoma = 010, Virginia = 001.

If the three political leaning values were one-hot encoded, they would be conservative = 100, moderate = 010, liberal = 001, but drop-first drops the first digit, giving conservative = 00, moderate = 10, liberal = 01.

If you use one-hot encoding, you include redundant information in a sense, in part because both one-hot encoding and drop-first encoding unambiguously identify a variable’s value, but drop-first does so using one less bit of information.

Now, if you train a linear regression model using a pseudo-inverse (either left-pseudo inverse via normal equations, or relaxed Moore-Penrose pseudo-inverse via SVD or QR decomposition), one-hot encoding creates multicollinearity in the data and the inverse operation will usually fail. So, when using a pseudo-inverse training technique, drop-first encoding is pretty much required.

If you train a linear regression model using stochastic gradient descent (SGD), you don’t run into the matrix inverse failure issue, and so one-hot encoding will work. However, because one-hot encoding has redundant information, SGD can produce wildly different results in the model weights and bias, but which give the same prediction accuracy. So, if you use one-hot encoding and SGD training, you get a model that predicts well, but one where the weights and bias can’t be interpreted well. Therefore, you might as well just use drop-first encoding.

Note that for neural networks, one-hot encoding is always used with SGD training because neural networks can’t be interpreted easily, so drop-first encoding provides no advantage.


It’s important to pay attention to details when working with machine learning. I love old science fiction movies from the 1950s. Minor characters can make a big difference in movies.

A relatively obscure actor named John Zaremba (1908-1986) played a supporting role in three of my favorites. He added a lot to the final result.

Left: In “Earth vs. the Flying Saucers” (1956), Zaremba plays Professor Kanter, a scientist who helped develop a sonic weapon to defeat an alien invasion.

Center: In “The Magnetic Monster” (1953), Zaremba plays Engineer Watson, the chief engineer for the Department of Power and Light. He receives the initial phone call alerting authorities to the existence of a new man-made element that doubles in mass every 11 hours, threatening to throw Earth out of orbit.

Right: In “20 Million Miles to Earth”, Zaremba plays Dr. Judson Uhl, a scientist involved with a secret manned mission to Venus. The spacecraft brings back a small Venusian creature who grows to giant size and causes havoc.


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

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

namespace LinearRegressionCategorical
{
  internal class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin C# linear regression ");
      Console.WriteLine("Predict salary from sex, age," +
        " State, politics ");

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

      string testFile =
        "C:\\VSM\\LinearRegressionCategorical\\" +
        "Data\\people_test_df.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("F5").
          PadLeft(8));

      // 2. create and train model
      double lrnRate = 0.001;
      int maxEpochs = 3000;
      int seed = 0;
      Console.WriteLine("\nSetting lrnRate = " +
        lrnRate.ToString("F4"));
      Console.WriteLine("Setting maxEpohcs = " +
        maxEpochs);

      Console.WriteLine("\nCreating and training" +
        " Linear Regression model using SGD ");
      LinearRegressor model =
        new LinearRegressor(seed);
      model.Train(trainX, trainY, lrnRate, maxEpochs);
      Console.WriteLine("Done ");

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

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

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

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

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

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

  } // class Program

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

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

    public void Train(double[][] trainX, double[] trainY,
      double lrnRate, int maxEpochs)
    {
      int n = trainX.Length;
      int dim = trainX[0].Length;
      this.weights = new double[dim];
      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];
      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)
        {
          int idx = indices[i];
          double[] x = trainX[idx];
          double predY = this.Predict(x);
          double actualY = trainY[idx];
          for (int j = 0; j "lt" dim; ++j)
            this.weights[j] -= lrnRate *
              (predY - actualY) * x[j];
          this.bias -= lrnRate * (predY - actualY);
        }
        if (epoch % (int)(maxEpochs / 5) == 0)
        {
          double mse = this.MSE(trainX, trainY);
          string s1 = "epoch = " + 
            epoch.ToString().PadLeft(5);
          string s2 = "  MSE = " + 
            mse.ToString("F6").PadLeft(8);
          Console.WriteLine(s1 + s2);
        }
      }
    }

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

    private static void Shuffle(int[] indices, Random rnd)
    {
      // Fisher-Yates
      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:

# people_train_df.txt
# sex (0 = male, 1 = female)
# age, state (michigan, nebraska, oklahoma), income,
# politics type (conservative, moderate, liberal)
#
1, 0.24, 0, 0, 0.29500, 0, 1
0, 0.39, 0, 1, 0.51200, 1, 0
1, 0.63, 1, 0, 0.75800, 0, 0
0, 0.36, 0, 0, 0.44500, 1, 0
1, 0.27, 1, 0, 0.28600, 0, 1
1, 0.50, 1, 0, 0.56500, 1, 0
1, 0.50, 0, 1, 0.55000, 1, 0
0, 0.19, 0, 1, 0.32700, 0, 0
1, 0.22, 1, 0, 0.27700, 1, 0
0, 0.39, 0, 1, 0.47100, 0, 1
1, 0.34, 0, 0, 0.39400, 1, 0
0, 0.22, 0, 0, 0.33500, 0, 0
1, 0.35, 0, 1, 0.35200, 0, 1
0, 0.33, 1, 0, 0.46400, 1, 0
1, 0.45, 1, 0, 0.54100, 1, 0
1, 0.42, 1, 0, 0.50700, 1, 0
0, 0.33, 1, 0, 0.46800, 1, 0
1, 0.25, 0, 1, 0.30000, 1, 0
0, 0.31, 1, 0, 0.46400, 0, 0
1, 0.27, 0, 0, 0.32500, 0, 1
1, 0.48, 0, 0, 0.54000, 1, 0
0, 0.64, 1, 0, 0.71300, 0, 1
1, 0.61, 1, 0, 0.72400, 0, 0
1, 0.54, 0, 1, 0.61000, 0, 0
1, 0.29, 0, 0, 0.36300, 0, 0
1, 0.50, 0, 1, 0.55000, 1, 0
1, 0.55, 0, 1, 0.62500, 0, 0
1, 0.40, 0, 0, 0.52400, 0, 0
1, 0.22, 0, 0, 0.23600, 0, 1
1, 0.68, 1, 0, 0.78400, 0, 0
0, 0.60, 0, 0, 0.71700, 0, 1
0, 0.34, 0, 1, 0.46500, 1, 0
0, 0.25, 0, 1, 0.37100, 0, 0
0, 0.31, 1, 0, 0.48900, 1, 0
1, 0.43, 0, 1, 0.48000, 1, 0
1, 0.58, 1, 0, 0.65400, 0, 1
0, 0.55, 1, 0, 0.60700, 0, 1
0, 0.43, 1, 0, 0.51100, 1, 0
0, 0.43, 0, 1, 0.53200, 1, 0
0, 0.21, 0, 0, 0.37200, 0, 0
1, 0.55, 0, 1, 0.64600, 0, 0
1, 0.64, 1, 0, 0.74800, 0, 0
0, 0.41, 0, 0, 0.58800, 1, 0
1, 0.64, 0, 1, 0.72700, 0, 0
0, 0.56, 0, 1, 0.66600, 0, 1
1, 0.31, 0, 1, 0.36000, 1, 0
0, 0.65, 0, 1, 0.70100, 0, 1
1, 0.55, 0, 1, 0.64300, 0, 0
0, 0.25, 0, 0, 0.40300, 0, 0
1, 0.46, 0, 1, 0.51000, 1, 0
0, 0.36, 0, 0, 0.53500, 0, 0
1, 0.52, 1, 0, 0.58100, 1, 0
1, 0.61, 0, 1, 0.67900, 0, 0
1, 0.57, 0, 1, 0.65700, 0, 0
0, 0.46, 1, 0, 0.52600, 1, 0
0, 0.62, 0, 0, 0.66800, 0, 1
1, 0.55, 0, 1, 0.62700, 0, 0
0, 0.22, 0, 1, 0.27700, 1, 0
0, 0.50, 0, 0, 0.62900, 0, 0
0, 0.32, 1, 0, 0.41800, 1, 0
0, 0.21, 0, 1, 0.35600, 0, 0
1, 0.44, 1, 0, 0.52000, 1, 0
1, 0.46, 1, 0, 0.51700, 1, 0
1, 0.62, 1, 0, 0.69700, 0, 0
1, 0.57, 1, 0, 0.66400, 0, 0
0, 0.67, 0, 1, 0.75800, 0, 1
1, 0.29, 0, 0, 0.34300, 0, 1
1, 0.53, 0, 0, 0.60100, 0, 0
0, 0.44, 0, 0, 0.54800, 1, 0
1, 0.46, 1, 0, 0.52300, 1, 0
0, 0.20, 1, 0, 0.30100, 1, 0
0, 0.38, 0, 0, 0.53500, 1, 0
1, 0.50, 1, 0, 0.58600, 1, 0
1, 0.33, 1, 0, 0.42500, 1, 0
0, 0.33, 1, 0, 0.39300, 1, 0
1, 0.26, 1, 0, 0.40400, 0, 0
1, 0.58, 0, 0, 0.70700, 0, 0
1, 0.43, 0, 1, 0.48000, 1, 0
0, 0.46, 0, 0, 0.64400, 0, 0
1, 0.60, 0, 0, 0.71700, 0, 0
0, 0.42, 0, 0, 0.48900, 1, 0
0, 0.56, 0, 1, 0.56400, 0, 1
0, 0.62, 1, 0, 0.66300, 0, 1
0, 0.50, 0, 0, 0.64800, 1, 0
1, 0.47, 0, 1, 0.52000, 1, 0
0, 0.67, 1, 0, 0.80400, 0, 1
0, 0.40, 0, 1, 0.50400, 1, 0
1, 0.42, 1, 0, 0.48400, 1, 0
1, 0.64, 0, 0, 0.72000, 0, 0
0, 0.47, 0, 0, 0.58700, 0, 1
1, 0.45, 1, 0, 0.52800, 1, 0
0, 0.25, 0, 1, 0.40900, 0, 0
1, 0.38, 0, 0, 0.48400, 0, 0
1, 0.55, 0, 1, 0.60000, 1, 0
0, 0.44, 0, 0, 0.60600, 1, 0
1, 0.33, 0, 0, 0.41000, 1, 0
1, 0.34, 0, 1, 0.39000, 1, 0
1, 0.27, 1, 0, 0.33700, 0, 1
1, 0.32, 1, 0, 0.40700, 1, 0
1, 0.42, 0, 1, 0.47000, 1, 0
0, 0.24, 0, 1, 0.40300, 0, 0
1, 0.42, 1, 0, 0.50300, 1, 0
1, 0.25, 0, 1, 0.28000, 0, 1
1, 0.51, 1, 0, 0.58000, 1, 0
0, 0.55, 1, 0, 0.63500, 0, 1
1, 0.44, 0, 0, 0.47800, 0, 1
0, 0.18, 0, 0, 0.39800, 0, 0
0, 0.67, 1, 0, 0.71600, 0, 1
1, 0.45, 0, 1, 0.50000, 1, 0
1, 0.48, 0, 0, 0.55800, 1, 0
0, 0.25, 1, 0, 0.39000, 1, 0
0, 0.67, 0, 0, 0.78300, 1, 0
1, 0.37, 0, 1, 0.42000, 1, 0
0, 0.32, 0, 0, 0.42700, 1, 0
1, 0.48, 0, 0, 0.57000, 1, 0
0, 0.66, 0, 1, 0.75000, 0, 1
1, 0.61, 0, 0, 0.70000, 0, 0
0, 0.58, 0, 1, 0.68900, 1, 0
1, 0.19, 0, 0, 0.24000, 0, 1
1, 0.38, 0, 1, 0.43000, 1, 0
0, 0.27, 0, 0, 0.36400, 1, 0
1, 0.42, 0, 0, 0.48000, 1, 0
1, 0.60, 0, 0, 0.71300, 0, 0
0, 0.27, 0, 1, 0.34800, 0, 0
1, 0.29, 1, 0, 0.37100, 0, 0
0, 0.43, 0, 0, 0.56700, 1, 0
1, 0.48, 0, 0, 0.56700, 1, 0
1, 0.27, 0, 1, 0.29400, 0, 1
0, 0.44, 0, 0, 0.55200, 0, 0
1, 0.23, 1, 0, 0.26300, 0, 1
0, 0.36, 1, 0, 0.53000, 0, 1
1, 0.64, 0, 1, 0.72500, 0, 0
1, 0.29, 0, 1, 0.30000, 0, 1
0, 0.33, 0, 0, 0.49300, 1, 0
0, 0.66, 1, 0, 0.75000, 0, 1
0, 0.21, 0, 1, 0.34300, 0, 0
1, 0.27, 0, 0, 0.32700, 0, 1
1, 0.29, 0, 0, 0.31800, 0, 1
0, 0.31, 0, 0, 0.48600, 1, 0
1, 0.36, 0, 1, 0.41000, 1, 0
1, 0.49, 1, 0, 0.55700, 1, 0
0, 0.28, 0, 0, 0.38400, 0, 0
0, 0.43, 0, 1, 0.56600, 1, 0
0, 0.46, 1, 0, 0.58800, 1, 0
1, 0.57, 0, 0, 0.69800, 0, 0
0, 0.52, 0, 1, 0.59400, 1, 0
0, 0.31, 0, 1, 0.43500, 1, 0
0, 0.55, 0, 0, 0.62000, 0, 1
1, 0.50, 0, 0, 0.56400, 1, 0
1, 0.48, 1, 0, 0.55900, 1, 0
0, 0.22, 0, 1, 0.34500, 0, 0
1, 0.59, 0, 1, 0.66700, 0, 0
1, 0.34, 0, 0, 0.42800, 0, 1
0, 0.64, 0, 0, 0.77200, 0, 1
1, 0.29, 0, 1, 0.33500, 0, 1
0, 0.34, 1, 0, 0.43200, 1, 0
0, 0.61, 0, 0, 0.75000, 0, 1
1, 0.64, 0, 1, 0.71100, 0, 0
0, 0.29, 0, 0, 0.41300, 0, 0
1, 0.63, 1, 0, 0.70600, 0, 0
0, 0.29, 1, 0, 0.40000, 0, 0
0, 0.51, 0, 0, 0.62700, 1, 0
0, 0.24, 0, 1, 0.37700, 0, 0
1, 0.48, 1, 0, 0.57500, 1, 0
1, 0.18, 0, 0, 0.27400, 0, 0
1, 0.18, 0, 0, 0.20300, 0, 1
1, 0.33, 1, 0, 0.38200, 0, 1
0, 0.20, 0, 1, 0.34800, 0, 0
1, 0.29, 0, 1, 0.33000, 0, 1
0, 0.44, 0, 1, 0.63000, 0, 0
0, 0.65, 0, 1, 0.81800, 0, 0
0, 0.56, 0, 0, 0.63700, 0, 1
0, 0.52, 0, 1, 0.58400, 1, 0
0, 0.29, 1, 0, 0.48600, 0, 0
0, 0.47, 1, 0, 0.58900, 1, 0
1, 0.68, 0, 0, 0.72600, 0, 1
1, 0.31, 0, 1, 0.36000, 1, 0
1, 0.61, 1, 0, 0.62500, 0, 1
1, 0.19, 1, 0, 0.21500, 0, 1
1, 0.38, 0, 1, 0.43000, 1, 0
0, 0.26, 0, 0, 0.42300, 0, 0
1, 0.61, 1, 0, 0.67400, 0, 0
1, 0.40, 0, 0, 0.46500, 1, 0
0, 0.49, 0, 0, 0.65200, 1, 0
1, 0.56, 0, 0, 0.67500, 0, 0
0, 0.48, 1, 0, 0.66000, 1, 0
1, 0.52, 0, 0, 0.56300, 0, 1
0, 0.18, 0, 0, 0.29800, 0, 0
0, 0.56, 0, 1, 0.59300, 0, 1
0, 0.52, 1, 0, 0.64400, 1, 0
0, 0.18, 1, 0, 0.28600, 1, 0
0, 0.58, 0, 0, 0.66200, 0, 1
0, 0.39, 1, 0, 0.55100, 1, 0
0, 0.46, 0, 0, 0.62900, 1, 0
0, 0.40, 1, 0, 0.46200, 1, 0
0, 0.60, 0, 0, 0.72700, 0, 1
1, 0.36, 1, 0, 0.40700, 0, 1
1, 0.44, 0, 0, 0.52300, 1, 0
1, 0.28, 0, 0, 0.31300, 0, 1
1, 0.54, 0, 1, 0.62600, 0, 0

Test data:

# people_test_df.txt
#
0, 0.51, 0, 0, 0.61200, 1, 0
0, 0.32, 1, 0, 0.46100, 1, 0
1, 0.55, 0, 0, 0.62700, 0, 0
1, 0.25, 0, 1, 0.26200, 0, 1
1, 0.33, 0, 1, 0.37300, 0, 1
0, 0.29, 1, 0, 0.46200, 0, 0
1, 0.65, 0, 0, 0.72700, 0, 0
0, 0.43, 1, 0, 0.51400, 1, 0
0, 0.54, 1, 0, 0.64800, 0, 1
1, 0.61, 1, 0, 0.72700, 0, 0
1, 0.52, 1, 0, 0.63600, 0, 0
1, 0.30, 1, 0, 0.33500, 0, 1
1, 0.29, 0, 0, 0.31400, 0, 1
0, 0.47, 0, 1, 0.59400, 1, 0
1, 0.39, 1, 0, 0.47800, 1, 0
1, 0.47, 0, 1, 0.52000, 1, 0
0, 0.49, 0, 0, 0.58600, 1, 0
0, 0.63, 0, 1, 0.67400, 0, 1
0, 0.30, 0, 0, 0.39200, 0, 0
0, 0.61, 0, 1, 0.69600, 0, 1
0, 0.47, 0, 1, 0.58700, 1, 0
1, 0.30, 0, 1, 0.34500, 0, 1
0, 0.51, 0, 1, 0.58000, 1, 0
0, 0.24, 0, 0, 0.38800, 1, 0
0, 0.49, 0, 0, 0.64500, 1, 0
1, 0.66, 0, 1, 0.74500, 0, 0
0, 0.65, 0, 0, 0.76900, 0, 0
0, 0.46, 1, 0, 0.58000, 0, 0
0, 0.45, 0, 1, 0.51800, 1, 0
0, 0.47, 0, 0, 0.63600, 0, 0
0, 0.29, 0, 0, 0.44800, 0, 0
0, 0.57, 0, 1, 0.69300, 0, 1
0, 0.20, 0, 0, 0.28700, 0, 1
0, 0.35, 0, 0, 0.43400, 1, 0
0, 0.61, 0, 1, 0.67000, 0, 1
0, 0.31, 0, 1, 0.37300, 1, 0
1, 0.18, 0, 0, 0.20800, 0, 1
1, 0.26, 0, 1, 0.29200, 0, 1
0, 0.28, 0, 0, 0.36400, 0, 1
0, 0.59, 0, 1, 0.69400, 0, 1
Posted in Machine Learning | Leave a comment