James D. McCaffrey

I write code almost every day. Writing code is a skill that must be practiced, and I enjoy writing code. One morning before work, I figured I’d run the well-known Diabetes Dataset through a nearest neighbors regression model. I didn’t expect very good results, and that’s what happened.

But the nearest neighbors regression results are useful as a baseline for comparison with more powerful regression techniques: quadratic regression, kernel ridge regression, neural network regression, decision tree regression, bagging tree regression, random forest regression, adaptive boosting regression, and gradient boosting regression.

The 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 but I suspect male = 1, female = 2 because there are 235 1 values and 206 2 values).

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 diabetes 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
. . .

I split the 442-items into a 342-item training set and a 100-item test set.

Nearest neighbors regression is arguably the simplest of all regression techniques. Suppose you want to predict the y value for an input vector x of predictor values. You scan through a set of training data (sometimes called the reference data) and kind the k-closest data items, where the value of k is typically about 3 or 5. Then you fetch the associated k target values from the training data, compute their average, and that average is the predicted y value.

I implemented a nearest neighbors regression model using C#. The output of the demo is:

Begin k-NN regression using C#

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

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

First three train y:
  0.1510
  0.0750
  0.1410

Creating k-NN regression model, k=5
Done

Accuracy train (within 0.10): 0.2602
Accuracy test (within 0.10): 0.1700

MSE train = 0.0027
MSE test = 0.0041

Predicting for trainX[0]
Predicted y = 0.1828

Explaining for train[0]
X = [  0]   0.59  . . .  0.09 | y = 0.1510 | dist = 0.0000 | wt = 0.2000
X = [215]   0.56  . . .  0.12 | y = 0.2630 | dist = 0.0598 | wt = 0.2000
X = [271]   0.59  . . .  0.09 | y = 0.1270 | dist = 0.0663 | wt = 0.2000
X = [341]   0.57  . . .  0.09 | y = 0.2630 | dist = 0.0678 | wt = 0.2000
X = [  8]   0.60  . . .  0.09 | y = 0.1100 | dist = 0.0686 | wt = 0.2000

Predicted y = 0.1828

End demo

The demo shows that nearest neighbors regression is highly interpretable. With k set to 5, to predict the y value for the first training item, the model finds the 5 closest training items: items [0], [215], [271], [341], [8]. The associated y values are (0.1510, 0.2630, 0.1270, 0.2630, 0.1100). The predicted y value is the average: (0.1510 + . . + 0.1100) / 5 = 0.1828. The demo weights each of the 5 items equally (wt = 0.2000). An alternative design is to give higher weights to closer items.

I didn’t expect great results and that’s what the model produced. The model scored just 26.02% accuracy on the training data (89 out of 342 correct) and 17.00% on the test data (17 out of 100 correct). A prediction is scored correct if it is within 10% of the true target value.

Interesting. It’s usually impossible to tell in advance if a nearest neighbors regression model will fit a particular dataset, and this is a good example of that.

Quadratic regression is an old machine learning technique that can often be effective (in terms of predictive accuracy). Shadow photography is an old technique that can often be effective (in terms of artistry, at least according to my eye).

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 NearestNeighborsRegression
{
  internal class NearestNeighborsProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin k-NN regression" +
        " using C# ");

      // 1. load data
      Console.WriteLine("\nLoading diabetes train" +
        " (342) and test (100) 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/fit model
      int numNeighbors = 5;
      Console.WriteLine("\nCreating k-NN regression " +
        "model, k=" + numNeighbors);
      KNNR model = new KNNR(numNeighbors);
      model.Train(trainX, trainY); // aka Fit()
      Console.WriteLine("Done ");

      // 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
      Console.WriteLine("\nPredicting for trainX[0] ");
      double[] x = trainX[0];
      double y = model.Predict(x);
      Console.WriteLine("Predicted y = " +
        y.ToString("F4").PadLeft(4));

      Console.WriteLine("\nExplaining for train[0] ");
      model.Explain(x);

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

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

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

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

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

  } // class Program

  public class KNNR
  {
    public int nNeighbors;
    public double[] weights;
    public double[][] trainX;
    public double[] trainY;

    public KNNR(int nNeighbors)
    {
      this.nNeighbors = nNeighbors;
      this.weights = new double[nNeighbors];
      for (int k = 0; k "lt" nNeighbors; ++k)
        this.weights[k] = 1.0 / nNeighbors;
      this.trainX = new double[0][]; // quasi-null
      this.trainY = new double[0];
    }

    public void Train(double[][] trainX,
      double[] trainY)
    {
      this.trainX = trainX; // store by ref
      this.trainY = trainY;
    }

    public double Predict(double[] x)
    {
      // compute distances from x to all trainX
      int N = trainX.Length;
      double[] dists = new double[N];
      for (int i = 0; i "lt" N; ++i)
        dists[i] = Distance(x, this.trainX[i]);
      // determine sort order
      int[] sortedIdxs = ArgSort(dists);
      double sum = 0.0;
      for (int k = 0; k "lt" this.nNeighbors; ++k)
      {
        int idx = sortedIdxs[k]; // near to far
        sum += this.trainY[idx] * this.weights[k];
      }
      return sum;
    }

    public void Explain(double[] x)
    {
      // 0. set up ordering/indices
      int n = this.trainX.Length;
      int[] indices = new int[n];
      for (int i = 0; i "lt" n; ++i)
        indices[i] = i;

      // 1. compute distances from x to all trainX
      double[] distances = new double[n];
      for (int i = 0; i "lt" n; ++i)
        distances[i] = Distance(x, this.trainX[i]);

      // 2. sort distances indices of X and Y, by distances
      Array.Sort(distances, indices);

      // 3. compute predicted y
      double predY = 0.0;
      for (int k = 0; k "lt" this.nNeighbors; ++k)
        predY += this.weights[k] * this.trainY[indices[k]];

      // 4. show info 
      for (int k = 0; k "lt" this.nNeighbors; ++k)
      {
        int i = indices[k];
        Console.Write("X = ");
        Console.Write("[" + i.ToString().PadLeft(3) + "] ");
 
        for (int j = 0; j "lt" this.trainX[i].Length; ++j)
          Console.Write(this.trainX[i][j].
            ToString("F2").PadLeft(6));

        Console.Write(" | y = ");
        Console.Write(this.trainY[i].ToString("F4"));
        Console.Write(" | dist = ");
        Console.Write(distances[k].ToString("F4"));
        Console.Write(" | wt = ");
        Console.Write(this.weights[k].ToString("F4"));
        Console.WriteLine("");
      }

      Console.WriteLine("\nPredicted y = " +
        predY.ToString("F4"));
    } // Explain

    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 actualY = dataY[i];
        if (Math.Abs(predY - actualY) "lt" 
          Math.Abs(pctClose * actualY))
          ++nCorrect;
        else
          ++nWrong;
      }
      return (nCorrect * 1.0) / (nCorrect + nWrong);
    }

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

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

    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);
      return indices;
    }
  } // class KNNR

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

One year ago, I wrote an article for the now-defunct Microsoft MSDN Magazine titled “Roach Infestation Optimization”. The C# code associated with that article vanished into the void, so I figured I’d rehydrate the demo program.

Roach infestation optimization (RIO) is a numerical optimization algorithm that’s loosely based on the behavior of common cockroaches such as Periplaneta americana. Let me explain.

In machine learning, a numerical optimization algorithm is often used to find a set of values for variables (usually called weights) that minimize some measure of error. The process of determining the values of the weights is called training the model. The idea is to use a collection of training data that has known correct output values. A numerical optimization algorithm is used to find the values of the weights that minimize the error between computed output values and known correct output values.

The most common training algorithms are based on calculus derivatives (stochastic gradient descent), but there are also algorithms that are based on the behaviors of natural systems. These are sometimes called bio-inspired algorithms. These bio-inspired algorithms are usually not practical, but they’re interesting.

The goal of the demo program is to use RIO to find the minimum value of the Rastrigin function with eight input variables. The Rastrigin function is a standard benchmark function used to evaluate the effectiveness of numerical optimization algorithms. The function has a known minimum value of 0.0 located at x = (0, 0, . . 0) where the number of 0 values is equal to the number of input values.

The Rastrigin function is extremely difficult for most optimization algorithms because it has many peaks and valleys that create local minimum values that can trap the algorithms. It’s not possible to easily visualize the Rastrigin function with eight input values, but you can get a good idea of the function’s characteristics by examining a graph of the function for two input values:

The output of the demo is:

Begin roach infestation optimization demo

Goal is minimize Rastrigin's function in 8 dimensions
Problem known min value = 0.0 at (0, 0, .. 0)

Setting number of roaches  = 50
Setting maximum time = 10000

Starting roach optimization

time =    500 best error = 104.76530
time =   1000 best error = 67.36795
time =   1500 best error = 12.11031
time =   2000 best error = 8.90949
time =   2500 best error = 3.99195
Mass extinction at t =   2500
time =   3000 best error = 3.99195
time =   3500 best error = 1.80840
time =   4000 best error = 1.02166
time =   4500 best error = 0.99382
time =   5000 best error = 0.00000
Mass extinction at t =   5000
time =   5500 best error = 0.00000
time =   6000 best error = 0.00000
time =   6500 best error = 0.00000
time =   7000 best error = 0.00000
time =   7500 best error = 0.00000
Mass extinction at t =   7500
time =   8000 best error = 0.00000
time =   8500 best error = 0.00000
time =   9000 best error = 0.00000
time =   9500 best error = 0.00000

Roach algorithm completed

Best error found = 0.000000 at:

0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000

End roach infestation optimization demo

The demo program sets the number of roaches to 50. Each simulated roach has a position that represents a possible solution to the minimization problem. More roaches increase the chance of finding the true optimal solution at the expense of performance.

RIO is an iterative process and requires a maximum loop counter value. The demo sets the maximum value to 10,000 iterations. The maximum number of iterations will vary from problem to problem.

In the demo, the best (smallest) error associated with the best roach position found so far was displayed every 500 time units. After the algorithm finished, the best position found for any roach was x = (0, 0, 0, 0, 0, 0, 0, 0), which is, in fact, the correct answer. But notice if the maximum number of iterations had been set to 3,000 instead of 10,000, RIO would not have found the one global minimum. RIO is not guaranteed to find an optimal solution in practical scenarios.

RIO is essentially a variation of the similar particle swarm optimization (PSO). The main difference is that in RIO, roaches that are close to each other will sometimes exchange information with each other.

Bio-inspired algorithms such as roach infestation optimization and particle swarm optimization are usually too slow to be practical. But if/when quantum computing becomes a reality, these algorithms could become essential techniques.

“Terra Formars” (2016) is an absolutely bonkers Japanese science fiction movie. In the 21st century, Earth plans to colonize Mars. To terraform the planet, the surface of Mars is seeded with moss and cockroaches to spread the moss.

This movie is one of many things Japanese that absolutely baffle me.

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

using System;

// "Roach Infestation Optimization" (2008)
// IEEE Swarm Intelligence Symposium
// T. Havens, C. Spain, N. Salmon, J. Keller

namespace RoachOptimization
{
  class RoachOptimizationProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin roach infestation " +
        "optimization demo\n");

      int dim = 8;
      int numRoaches = 50;
      int tMax = 10000;
      Random rnd = new Random(1);
      
      Console.WriteLine("Goal is minimize Rastrigin's " +
        "function in " + dim + " dimensions");
      Console.WriteLine("Problem known min value = 0.0 " +
        "at (0, 0, .. 0) \n");
      Console.WriteLine("Setting number of roaches  = " +
        numRoaches);
      Console.WriteLine("Setting maximum time = " +
        tMax);
      
      Console.WriteLine("\nStarting roach optimization \n");
      double[] answer = 
        SolveRastrigin(dim, numRoaches, tMax, rnd);
      Console.WriteLine("\nRoach algorithm completed ");

      double err = Error(answer);
      Console.WriteLine("\nBest error found = " +
        err.ToString("F6") + " at: \n");

      for (int i = 0; i "lt" dim; ++i)
        Console.Write(answer[i].ToString("F4") + " ");
      Console.WriteLine("");

      Console.WriteLine("\nEnd roach infestation " +
        "optimization demo ");
      Console.ReadLine();
    } // Main

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

    static double[] SolveRastrigin(int dim, int numRoaches,
      int tMax, Random rnd)
    {
      // estimate solution to Rastrigin's function
      // using roach infestation optimization
      double C0 = 0.7;  // RIO is based on PSO
      double C1 = 1.43;
      // probs of exchanging info based on number neighbors
      // actual roach behavior: 0.49, 0.63, 0.65
      // double[] A = new double[] { 0.49, 0.63, 0.65 }; 
      double[] A = new double[] { 0.2, 0.3, 0.4 }; // better 

      // smaller divisor (e.g., 5 instead of 10)
      // roaches stay longer
      int tHunger = tMax / 10;  
      double minX = -10.0;
      double maxX = 10.0;
      int extinction = tMax / 4;  // mass extinction restart
      //Random rnd = new Random(seed);  // random order

      Roach[] herd = new Roach[numRoaches];  // 'intrusion'
      for (int i = 0; i "lt" numRoaches; ++i)
        herd[i] = new Roach(dim, minX, maxX, tHunger, i);

      int t = 0;  // loop counter (time)
      int[] indices = new int[numRoaches];
      for (int i = 0; i "lt" numRoaches; ++i)
        indices[i] = i;

      double bestError = double.MaxValue;  // by any roach
      double[] bestPosition = new double[numRoaches];

      int displayInterval = tMax / 20; // 20 times

      // distances between all pairs of roaches
      double[][] distances = new double[numRoaches][];
      for (int i = 0; i "lt" numRoaches; ++i)
        distances[i] = new double[numRoaches];

      while (t "lt" tMax) // main loop
      {
        if (t "gt" 0 && t % displayInterval == 0)
        {
          Console.Write("time = " +
            t.ToString().PadLeft(6));
          Console.WriteLine(" best error = " +
            bestError.ToString("F5"));
        }

        // compute distances between all pairs
        for (int i = 0; i "lt" numRoaches - 1; ++i)
        {
          for (int j = i + 1; j "lt" numRoaches; ++j)
          {
            double d = 
              Distance(herd[i].position, herd[j].position);
            distances[i][j] = distances[j][i] = d;
          }
        }

        // find 'quarter' median distance
        double[] sortedDists = 
          new double[numRoaches * (numRoaches - 1) / 2];
        int k = 0;
        for (int i = 0; i "lt" numRoaches - 1; ++i)
        {
          for (int j = i + 1; j "lt" numRoaches; ++j)
          {
            sortedDists[k++] = distances[i][j];
          }
        }

        Array.Sort(sortedDists);
        double medianDist = 
          sortedDists[(int)(sortedDists.Length / 4)];

        Shuffle(indices, rnd); // process roaches rnd order

        for (int i = 0; i "lt" numRoaches; ++i)  // each roach
        {
          int idx = indices[i]; // roach index

          Roach curr = herd[idx]; // ref to current roach
          int numNeighbors = 0;

          // calculate number of neighbors of curr roach
          for (int j = 0; j "lt" numRoaches; ++j)  
          {
            if (j == idx) continue; // not a neighbor self
            double d = distances[idx][j];
            if (d "lt" medianDist) // is a neighbor to roach[j]
              ++numNeighbors;
          }

          // at most 3 neighbors are relevant
          int effectiveNeighbors = numNeighbors;
          if (effectiveNeighbors "gte" 3)  
            effectiveNeighbors = 3;

          // possibly swap group best information
          for (int j = 0; j "lt" numRoaches; ++j)
          {
            if (j == idx) continue; // not neighbor of self
            if (effectiveNeighbors == 0) continue;
            double prob = rnd.NextDouble();
            if (prob "gt" A[effectiveNeighbors - 1])
              continue; // don't always swap

            double d = distances[idx][j];
            if (d "lt" medianDist) // is a neighbor
            {
              // exchange info if curr roach better
              if (curr.error "lt" herd[j].error) // 
              {
                for (int p = 0; p "lt" dim; ++p)
                {
                  // use curr roach's personal best
                  // as the group best for both roaches
                  herd[j].groupBestPos[p] = 
                    curr.personalBestPos[p];  
                  curr.groupBestPos[p] = 
                    curr.personalBestPos[p];
                }
              }
              else // roach[j] is better than curr
              {
                for (int p = 0; p "lt" dim; ++p)
                {
                  curr.groupBestPos[p] = 
                    herd[j].personalBestPos[p];
                  herd[j].groupBestPos[p] = 
                    herd[j].personalBestPos[p];
                }
              }

            } // if a neighbor
          } // j, each neighbor

          if (curr.hunger "lt" tHunger) // move roach
          {
            // new velocity
            for (int p = 0; p "lt" dim; ++p)
              curr.velocity[p] = 
                (C0 * curr.velocity[p]) +
                (C1 * rnd.NextDouble() *
                (curr.personalBestPos[p] - 
                curr.position[p])) +
                (C1 * rnd.NextDouble() *
                (curr.groupBestPos[p] - 
                curr.position[p]));

            // update position
            for (int p = 0; p "lt" dim; ++p)
              curr.position[p] = 
                curr.position[p] + curr.velocity[p];

            double e = Error(curr.position); // Rastrigin
            curr.error = e;

            // new personal best position?
            if (curr.error "lt" curr.personalBestErr)
            {
              curr.personalBestErr = curr.error;
              for (int p = 0; p "lt" dim; ++p)
                curr.personalBestPos[p] = curr.position[p];
            }

            // new global best position?
            if (curr.error "lt" bestError)
            {
              bestError = curr.error;
              for (int p = 0; p "lt" dim; ++p)
                bestPosition[p] = curr.position[p];
            }

            ++curr.hunger;
          }
          else // hungry. leave group to random position
          {
            herd[idx] = 
              new Roach(dim, minX, maxX, tHunger, t);
          }

        } // each roach

        // mass extinction?
        if (t "gt" 0 && t % extinction == 0)  
        {
          Console.WriteLine("Mass extinction at t = " +
            t.ToString().PadLeft(6));
          for (int i = 0; i "lt" numRoaches; ++i)
            herd[i] = new Roach(dim, minX, maxX, tHunger, i);
        }

        ++t;
      } // main while loop

      return bestPosition;

    } // SolveRastrigin

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

    public static double Error(double[] x)
    {
      // Rastrigin function has min value of 0.0
      double trueMin = 0.0;
      double rastrigin = 0.0;
      for (int i = 0; i "lt" x.Length; ++i)
      {
        double xi = x[i];
        rastrigin += (xi * xi) -
          (10 * Math.Cos(2 * Math.PI * xi)) + 10;
      }

      return (rastrigin - trueMin) * (rastrigin - trueMin);
    }

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

    static double Distance(double[] pos1, double[] pos2)
    {
      double sum = 0.0;
      for (int i = 0; i "lt" pos1.Length; ++i)
        sum += (pos1[i] - pos2[i]) * (pos1[i] - pos2[i]);
      return Math.Sqrt(sum);
    }

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

    static void Shuffle(int[] indices, Random rnd)
    {
      //Random rnd = new Random(seed);
      for (int i = 0; i "lt" indices.Length; ++i)
      {
        int r = rnd.Next(i, indices.Length);
        int tmp = indices[r];
        indices[r] = indices[i];
        indices[i] = tmp;
      }
    }

  } // Program

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

  public class Roach
  {
    public int dim;
    public double[] position;
    public double[] velocity;
    public double error; // at curr position

    public double[] personalBestPos; // best position found
    public double personalBestErr; // 

    public double[] groupBestPos; // by any roach
 
    public int hunger;
    private Random rnd;

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

    public Roach(int dim, double minX, double maxX,
      int tHunger, int rndSeed)
    {
      this.dim = dim;
      this.position = new double[dim];
      this.velocity = new double[dim];
      this.personalBestPos = new double[dim];
      this.groupBestPos = new double[dim];

      this.rnd = new Random(rndSeed);
      this.hunger = this.rnd.Next(0, tHunger);

      for (int i = 0; i "lt" dim; ++i)
      {
        this.position[i] = 
          (maxX - minX) * rnd.NextDouble() + minX;
        this.velocity[i] = 
          (maxX - minX) * rnd.NextDouble() + minX;
        this.personalBestPos[i] = this.position[i];
        this.groupBestPos[i] = this.position[i];
      }

      this.error = 
        RoachOptimizationProgram.Error(this.position);
      this.personalBestErr = this.error;
    }

  } // class Roach

} // ns

If you have a matrix A, the singular value decomposition (SVD) gives you a matrix U, a vector s, and a matrix Vh, so that A = U * S * Vh. The s is a vector like [2, 4, 6], and S is a square matrix with the s values on the diagonal, 0s off the diagonal.

SVD is used in dozens of important software algorithms.

SVD is an extremely complicated topic, and computing the SVD of a matrix is one of the most challenging problems in all of computer programming. Many researchers worked for many years on the problem, and came up with many different versions of SVD implementations that vary in accuracy, complexity, and robustness. The standard SVD implementation in LAPACK is over 10,000 lines of Fortran code.

I put together an implementation of SVD (Householder + QR version) using the C# language. Unfortunately, my implementation is not as stable as my older implementation that uses the Jacobi algorithm. The first part of my demo output is:

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

Computing SVD
Done

U =
   0.319267   0.295047   0.358160   0.465731   0.652040
   0.625260  -0.614481   0.189087   0.174205  -0.138090
  -0.129989   0.032917  -0.342862  -0.160369   0.572690
   0.284056   0.494504  -0.490845   0.494308  -0.381366
   0.486033   0.495867   0.263713  -0.648653  -0.103790
   0.416300  -0.209426  -0.638701  -0.248874   0.267559

s =
  15.967661  12.793149   6.297366   5.706879   2.478679

Vh =
   0.296544   0.035215   0.708796  -0.130799   0.625557
   0.217254   0.416871   0.261754   0.803400  -0.255056
  -0.745282   0.555722  -0.030757   0.039094   0.365040
  -0.187161  -0.609726  -0.196995   0.579569   0.467438
   0.523820   0.379983  -0.623971   0.003540   0.438033

Behind the scenes, I verified my implementation by feeding the same source matrix A to the NumPy library np.linalg.svd() function to make sure I got the same results. I also constructed the S matrix from the s vector, and then verified that A = U * S * Vh.

For the second part of my demo, I used a “wide” matrix that has more columns than rows – the SVD algorithm works differently for tall matrices than it does for wide matrices:

Source (wide) matrix:
   3.0   1.0   1.0   0.0   5.0
  -1.0   3.0   1.0  -2.0   4.0
   0.0   2.0   2.0   1.0  -3.0

Computing SVD
Done

U =
  -0.727149   0.195146  -0.658158
  -0.621912  -0.593195   0.511219
   0.290654  -0.781049  -0.552705

s =
   7.639218   4.023860   3.232785

Vh =
  -0.204149  -0.263322  -0.100502   0.200868  -0.915716
   0.292911  -0.781970  -0.487131   0.100734   0.235122
  -0.768902  -0.071118  -0.387390  -0.487240   0.127506

My C# implementation is intended mostly for machine learning scenarios where accuracy to 1.0e-8 is good enough, as opposed to scientific programming where accuracy to 1.0e-15 is typical.

I worked on this SVD implementation off-and-on for several years. I finally got some extended uninterrupted time and coded up a demo but the implementation took several weeks to write. I used over a dozen major resources. The implementation here has all normal error-checking removed to keep the number of lines of code as small as possible.

I’ve written some complex code over the years but this SVD implementation is one of the most complex. I consider the code presented in this blog post to be one of the interesting explorations of my software engineering career. Unfortunately, this implementation isn’t as stable as my older SVD using the Jacobi algorithm.

The development of an algorithm to compute the SVD of a matrix was one of the key advances in the history of computer programming.

Previous innovations included first electrically powered bumpers (1948), modern flipper geometry (1950), automatic reel scoring (1953), and drop-targets (1962).

Demo program. Very long, very complicated. 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;

// no significant .NET dependencies

namespace MatrixSingularValueDecomposition
{
  internal class SVDProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin SVD using scratch C# ");
      double[][] A = new double[6][];
      A[0] = new double[] { 1, 2, 3, 4, 5 };
      A[1] = new double[] { 0, -3, 5, -7, 9 };
      A[2] = new double[] { 2, 0, -2, 0, -2 };
      A[3] = new double[] { 4, -1, 5, 6, 1 };
      A[4] = new double[] { 3, 6, 8, 2, 2 };
      A[5] = new double[] { 5, -2, 4, -4, 3 };

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

      Console.WriteLine("\nComputing SVD ");
      double[][] U;
      double[] s;
      double[][] Vh;
      SVD.MatDecomp(A, out U, out s, out Vh);
      Console.WriteLine("Done ");

      Console.WriteLine("\nU = ");
      MatShow(U, 6, 11);

      Console.WriteLine("\ns = ");
      for (int i = 0; i "lt" s.Length; ++i)
        Console.Write(s[i].ToString("F6").PadLeft(11));
      Console.WriteLine("");

      Console.WriteLine("\nVh = ");
      MatShow(Vh, 6, 11);

      double[][] B = new double[3][];
      B[0] = new double[] { 3, 1, 1, 0, 5 };
      B[1] = new double[] { -1, 3, 1, -2, 4 };
      B[2] = new double[] { 0, 2, 2, 1, -3 };

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

      Console.WriteLine("\nSource (wide) matrix: ");
      MatShow(B, 1, 6);

      Console.WriteLine("\nComputing SVD ");
      SVD.MatDecomp(B, out U, out s, out Vh);
      Console.WriteLine("Done ");

      Console.WriteLine("\nU = ");
      MatShow(U, 6, 11);

      Console.WriteLine("\ns = ");
      for (int i = 0; i "lt" s.Length; ++i)
        Console.Write(s[i].ToString("F6").PadLeft(11));
      Console.WriteLine("");

      Console.WriteLine("\nVh = ");
      MatShow(Vh, 6, 11);

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

    // helpers for Main()

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

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

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

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      // ex: double[][] A =
      //  MatLoad("data.txt", new int[] {0,1,2}, ',', "#");
      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();
    }

  } // class Program

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

  public class SVD
  {
    public static void MatDecomp(double[][] A,
      out double[][] U, out double[] s, out double[][] Vt,
      double tol = 1.0e-12)
    {
      int m = A.Length; int n = A[0].Length;
      if (m "gte" n)
      {
        double[][] Q; double[][] R;
        MatDecompQR(A, out Q, out R);

        double[][] RtR = MatProduct(MatTranspose(R), R);

        double[] evals; double[][] V;
        MatDecompEigen(RtR, out evals, out V);

        int[] sortIdxs = VecReversed(ArgSort(evals));
        evals = VecSortedUsingIdxs(evals, sortIdxs);
        V = MatSortedUsingIdxs(V, sortIdxs);

        double[] ss = VecSqrt(VecClip(evals, 0.0));
        int[] idxs = VecIndicesWhereNonNegative(ss, tol);
        double[] ssNonZero = VecRemoveZeros(ss, idxs);
        double[][] VNonZero = MatRemoveCols(V, idxs);

        double[][] tmp1 = MatProduct(R, VNonZero);
        double[][] tmp2 = MatDivVector(tmp1, ssNonZero);
        double[][] UU = MatProduct(Q, tmp2);

        U = UU;
        s = ssNonZero;
        Vt = MatTranspose(VNonZero);
      } // m "gte" n
      else
      {
        // use the transpose trick
        double[][] dummyU;
        double[] dummyS;
        double[][] dummyVt;
        MatDecomp(MatTranspose(A),
          out dummyU, out dummyS, out dummyVt);
        U = MatTranspose(dummyVt);
        s = dummyS;
        Vt = MatTranspose(dummyU);
      } // m "lt" n

    } // MatDecomp()

    // ------------------------------------------------------
    // primary helpers: MatDecompQR() and MatDecompEigen()
    // ------------------------------------------------------

    private static void MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // reduced QR decomp using Householder reflections
      int m = A.Length; int n = A[0].Length;
      double[][] QQ = MatIdentity(m);  // working Q
      double[][] RR = MatCopy(A);  // working R

      for (int k = 0; k "lt" n; ++k)
      {
        double[] x = MatColToEnd(RR, k); // RR[k:, k]
        double normx = VecNorm(x);
        if (Math.Abs(normx) "lt" 1.0e-12) continue;

        double[] v = VecCopy(x);
        double sign;
        if (x[0] "gte" 0.0) sign = 1.0; else sign = -1.0;
        v[0] += sign * normx;
        double normv = VecNorm(v);
        for (int i = 0; i "lt" v.Length; ++i)
          v[i] /= normv;

        // apply reflection to R
        double[][] tmp1 = 
          MatRowColToEnd(RR, k); // RR[k:, k:]
        double[] tmp2 = VecMatProduct(v, tmp1);
        double[][] tmp3 = VecOuter(v, tmp2);
        double[][] B = MatScalarMult(2.0, tmp3);
        MatSubtractCorner(RR, k, B);  // R[k:, k:] -= B

        // accumulate Q
        double[][] tmp4 = 
          MatExtractAllRowsColToEnd(QQ, k); // QQ[:, k:]
        double[] tmp5 = MatVecProduct(tmp4, v);
        double[][] tmp6 = VecOuter(tmp5, v);
        double[][] C = MatScalarMult(2.0, tmp6);
        MatSubtractRows(QQ, k, C); // QQ[:, k:] -= C
      } // for k

      // return parts of QQ and RR
      Q = MatExtractFirstCols(QQ, n); // Q[:, :n]
      R = MatExtractFirstRows(RR, n); // R[:n, :]
    } // MatDecompQR()

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

    private static void MatDecompEigen(double[][] A,
      out double[] evals, out double[][] Evecs,
      double tol = 1.0e-12, int maxIter = 5000)
    {
      // eigen-decomp of symmetric matrix using QR
      int m = A.Length; int n = A[0].Length;
      double[][] AA = MatCopy(A);  // don't destroy A
      double[][] V = MatIdentity(n);

      int iter = 0;
      while (iter "lt" maxIter)
      {
        double[][] Q;
        double[][] R;
        MatDecompQR(AA, out Q, out R);
        double[][] Anew = MatProduct(R, Q);
        V = MatProduct(V, Q);

        // # check convergence (off-diagonal norm)
        double[] tmp1 = MatDiag(Anew);
        double[][] tmp2 = MatCreateDiag(tmp1);
        double[][] tmp3 = MatSubtract(Anew, tmp2);
        double norm = MatNorm(tmp3);
        if (norm "lt" tol)
        {
          AA = Anew; break;
        }
        AA = Anew;
        ++iter;
      } // iter

      if (iter == maxIter)
        Console.WriteLine("Warn: exceeded maxIter ");

      evals = MatDiag(AA);
      Evecs = V;
    } // MatDecompEigen()

    // ------------------------------------------------------
    // 32 secondary helpers
    // ------------------------------------------------------

    private static int[] ArgSort(double[] vec)
    {
      int n = vec.Length;
      int[] idxs = new int[n];
      for (int i = 0; i "lt" n; ++i)
        idxs[i] = i;

      double[] dup = new double[n];
      for (int i = 0; i "lt" n; ++i)
        dup[i] = vec[i]; // don't modify vec

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

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

    private static double[][] MatMake(int nr, int nc)
    {
      double[][] result = new double[nr][];
      for (int i = 0; i "lt" nr; ++i)
        result[i] = new double[nc];
      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[][] MatCopy(double[][] A)
    {
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[i][j];
      return result;
    }

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

    private static double[] MatDiag(double[][] A)
    {
      // return diag elements of A as a vector
      int m = A.Length; int n = A[0].Length; // m == n
      double[] result = new double[m];
      for (int i = 0; i "lt" m; ++i)
        result[i] = A[i][i];
      return result;
    }

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

    private static double[][] MatCreateDiag(double[] v)
    {
      // create a matrix using v as diagonal elements
      int n = v.Length;
      double[][] result = MatMake(n, n);
      for (int i = 0;i "lt" n; ++i)
        result[i][i] = v[i];
      return result;
    }

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

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

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

    private static double[][] MatSubtract(double[][] A,
      double[][] B)
    {
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[i][j] - B[i][j];
      return result;
    }

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

    private static double[] MatColToEnd(double[][] A,
      int k)
    {
      // last part of col k, from [k,k] to end
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[m - k];
      for (int i = 0; i "lt" m - k; ++i)
        result[i] = A[i + k][k];
      return result;
    }

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

    private static double[][] MatRowColToEnd(double[][] A,
      int k)
    {
      // block of A, row k to end, col k to end
      // A[k:, k:]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m-k, n-k);
      for (int i = 0; i "lt" m - k; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          result[i][j] = A[i + k][j + k];
      return result;
    }

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

    private static double[] MatVecProduct(double[][] A,
      double[] v)
    {
      // A * v
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[m];
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i] += A[i][j] * v[j];
      return result;
    }

    private static double[] VecMatProduct(double[] v,
      double[][] A)
    {
      // v * A
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[n];
      for (int j = 0; j "lt" n; ++j)
        for (int i = 0; i "lt" m; ++i)
          result[j] += v[i] * A[i][j];
      return result;
    }

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

    private static double[][] MatScalarMult(double x,
      double[][] A)
    {
      // x * A element-wise
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = x * A[i][j];
      return result;
    }

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

    private static void MatSubtractCorner(double[][] A,
      int k, double[][] C)
    {
      // subtract elements in C from lower right A
      // A[k:, k:] -= C
      int m = A.Length; int n = A[0].Length;
      for (int i = 0; i "lt" m - k; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          A[i + k][j + k] -= C[i][j];

      return;  // modify A in-place 
    }

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

    private static void MatSubtractRows(double[][] A,
      int k, double[][] C)
    {
      // A[:, k:] -= C
      int m = A.Length; int n = A[0].Length;
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n-k; ++j)
          A[i][j+k] -= C[i][j];

      return;  // modify A in-place
    }

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

    private static double[][] 
      MatExtractAllRowsColToEnd(double[][] A, int k)
    {
      // A[:, k:]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n - k);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          result[i][j] = A[i][j + k];
      return result;
    }

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

    private static double[][] 
      MatExtractFirstCols(double[][] A, int k)
    {
      // A[:, :n]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, k);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" k; ++j)
          result[i][j] = A[i][j];
      return result;
    }

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

    private static double[][] 
      MatExtractFirstRows(double[][] A, int k)
    {
      // A[:n, :]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(k, n);
      for (int i = 0; i "lt" k; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[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;
    }

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

    private static double[][] MatRemoveCols(double[][] A,
      int[] idxs)
    {
      // remove the columns specified in idxs
      int m = A.Length; int n = A[0].Length;
      int countNonzeos = 0;
      for (int i = 0; i "lt" idxs.Length; ++i)
        if (idxs[i] == 1) ++countNonzeos;

      double[][] result = MatMake(m, countNonzeos);
      int k = 0; // destination col
      for (int j = 0; j "lt" n; ++j)
      {
        if (idxs[j] == 0) continue; // skip this col
        for (int i = 0; i "lt" m; ++i)
        {
          result[i][k] = A[i][j];
        }
        ++k;
      }
      return result;
    }

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

    private static double[][] MatDivVector(double[][] A,
      double[] v)
    {
      // v must be all non-zero
      // divide each row of A by v, element-wise
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j= 0; j "lt" n; ++j)
          result[i][j] = A[i][j] / v[j];
      return result;
    }

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

    private static double MatNorm(double[][] A)
    {
      // treat A as one big vector
      int m = A.Length; int n = A[0].Length;
      double sum = 0.0;
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          sum += A[i][j] * A[i][j];
      return Math.Sqrt(sum); 
    }

    private static double VecNorm(double[] v)
    {
      // standard vector norm
      double sum = 0.0;
      for (int i = 0; i "lt" v.Length; ++i)
        sum += v[i] * v[i];
      return Math.Sqrt(sum);
    }

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

    private static double[][] 
      MatSortedUsingIdxs(double[][] A, int[] idxs)
    {
      // arrange columns of A using order in idxs
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          int srcCol = idxs[j];
          result[i][j] = A[i][srcCol];
        }

      }
      return result;
    }

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

    private static double[] 
      VecSortedUsingIdxs(double[] v, int[] idxs)
    {
      // arrange values in v using order specified in idxs
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = v[idxs[i]];
      return result;
    }

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

    private static double[] VecCopy(double[] v)
    {
      double[] result = new double[v.Length];
      for (int i = 0; i "lt" v.Length; ++i)
        result[i] = v[i];
      return result;
    }

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

    private static double[][] VecOuter(double[] v1,
      double[] v2)
    {
      // vector outer product
      double[][] result = MatMake(v1.Length, v2.Length);
      for (int i = 0; i "lt" v1.Length; ++i)
        for (int j = 0; j "lt" v2.Length; ++j)
          result[i][j] = v1[i] * v2[j];
      return result;
    }

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

    private static int[] VecReversed(int[] v)
    {
      // return vector of v elements in reversed order
      int n = v.Length;
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = v[n-1-i];
      return result;
    }

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

    private static double[] VecClip(double[] v, double minVal)
    {
      // used to make sure no negative values
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        if (v[i] "lt" minVal)
          result[i] = minVal;
        else
          result[i] = v[i];
      }
      return result;
    }

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

    private static double[] VecSqrt(double[] v)
    {
      // sqrt each element in v
      // assumes v has no negative values
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = Math.Sqrt(v[i]);
      return result;
    }

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

    private static int[] 
      VecIndicesWhereNonNegative(double[] v, double tol)
    {
      // indices of non-negative values (within tolerance)
      int n = v.Length;
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
      {
        if (Math.Abs(v[i]) "gt" tol)
          result[i] = 1;
      }
      return result;
    }

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

    private static double[] VecRemoveZeros(double[] v,
      int[] idxs)
    {
      // if idx = 1, val is non-zero, if 0 val is near zero
      int n = v.Length;
      List"lt"double"gt" lst = new List"lt"double"gt"();
      for (int i = 0; i "lt" n; ++i)
      {
        if (idxs[i] == 1)
          lst.Add(v[i]);
      }
      return lst.ToArray();
    }

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

  } // class SVD

} // ns

The condition number of a mathematical object is a number that is a measure of how sensitivity the object is to small changes. There are several types of condition numbers. In machine learning you might want to know how sensitive a square matrix A that defines a set of linear equations Ax = b is. Or you might want to know how sensitive a tall matrix (more rows than columns) that holds training data for a prediction model is.

The one-norm condition number of a matrix A is defined as cond(A) * cond(Ainv), in words, the one-norm of A times the one-norm of the inverse of A.

The one-norm of a matrix is the maximum of the column sums of the magnitudes of the values. For example, if A is:

   3.0   0.0  -2.0   5.0
  -1.0   4.0   6.0   3.0
   4.0   1.0   0.0   3.0
  -3.0   2.0   4.0   5.0

The column sums of the absolute values are 11, 7, 12, 16, and so the one-norm is 16.

Note: Another common condition number for matrices is the two-norm condition number. It is defines as the ratio of the largest singular value of the matrix divided by the smallest singular value. Computing singular values is one of the most difficult tasks in all of numerical programming.

I put together a demo using the C# language. To account for tall machine learning style matrices, I computed the pseudo-inverse (which applies to any shape matrix) instead of the standard inverse (which applies to only square matrices).

There are dozens of ways to compute a pseudo-inverse, all of them very complicated. I used QR decomposition via the modified Gram-Schmidt algorithm, which is, in my opinion, the simplest possible algorithm for pseudo-inverse.

The output of my demo:

Begin one-norm condition number demo

Source (square) matrix A:
   3.0   0.0  -2.0   5.0
  -1.0   4.0   6.0   3.0
   4.0   1.0   0.0   3.0
  -3.0   2.0   4.0   5.0

The one-norm condition number = 53.3333

Source (on-square) matrix B:
   3.0   0.0  -2.0   5.0
  -1.0   4.0   6.0   3.0
   4.0   1.0   0.0   3.0
  -3.0   2.0   4.0   5.0
   5.0   3.0  -1.0   0.0

The one-norm condition number = 16.5219

End demo

It’s very difficult to interpret condition numbers other than larger condition numbers mean more sensitive, in other words, small changes in the source matrix can lead to very large chances in the result (linear equation coefficients, matrix inverse for training, etc.)

I have a strange liking for comedy movies that have a special scene/condition: obviously fake bats. To me, these scenes are hilarious.

Right: In “Ace Ventura: When Nature Calls” (1995), detective Ace Ventura (actor Jim Carrey) loves all animals . . except bats. He is hired to go to Nibia, Africa to find a kidnapped sacred white bat. I’m not a Jim Carrey fan, but this movie has some hilarious scenes. My grade = B+.

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

// matrix 1-norm condition number

using System.IO;

namespace MatrixConditionNumberDefinition
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin one-norm condition" +
        " number demo ");

      // 1. set up source matrix
      double[][] A = new double[4][];
      A[0] = new double[] { 3, 0, -2, 5 };
      A[1] = new double[] { -1, 4, 6, 3 };
      A[2] = new double[] { 4, 1, 0, 3 };
      A[3] = new double[] { -3, 2, 4, 5 };
      // A[4] = new double[] { 5, 3, -1, 0 };

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

      double cnA = MatOneNormCondition(A);
      Console.WriteLine("\nThe one-norm condition " +
        "number = " + cnA.ToString("F4"));

      double[][] B = new double[5][];
      B[0] = new double[] { 3, 0, -2, 5 };
      B[1] = new double[] { -1, 4, 6, 3 };
      B[2] = new double[] { 4, 1, 0, 3 };
      B[3] = new double[] { -3, 2, 4, 5 };
      B[4] = new double[] { 5, 3, -1, 0 };

      Console.WriteLine("\nSource (on-square) matrix B: ");
      MatShow(B, 1, 6);

      double cnB = MatOneNormCondition(B);
      Console.WriteLine("\nThe one-norm condition " +
        "number = " + cnB.ToString("F4"));

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

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

    static double MatOneNormCondition(double[][] A)
    {
      // one-norm condition for a matrix in a
      // machine learning scenarios where nRows gte nCols
      // cond = norm(A) * norm(Ainv)
      double a = MatOneNorm(A);
      double b = MatOneNorm(QRGramSchmidt.MatPseudoInv(A));
      return a * b;
    }


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

    static double MatOneNorm(double[][] A)
    {
      // max col sum of abs values
      // do not assume A is a square matrix
      // assume m gte n
      int m = A.Length; int n = A[0].Length;
      double result = 0;
      for (int j = 0; j "lt" n; ++j) // each column
      {
        double colSum = 0.0;
        for (int i = 0; i "lt" m; ++i) // walk down column
        {
          colSum += Math.Abs(A[i][j]);
        }
        if (colSum "gt" result)
          result = colSum;
      }
      return result;
    }

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

    public static void MatShow(double[][] A, int dec,
      int wid)
    {
      int nRows = A.Length;
      int nCols = A[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 = A[i][j]; // avoid "-0.00000"
          if (Math.Abs(v) "lt" small) v = 0.0;
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

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

  } // class Program

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

  public class QRGramSchmidt
  {
    // container class for relaxed Moore-Penrose
    // pseudo-inverse using QR decomposition with modified
    // Gram-Schmidt algorithm

    public static double[][] MatPseudoInv(double[][] M)
    {
      // relaxed Moore-Penrose pseudo-inverse using QR-GS
      // A = Q*R, pinv(A) = inv(R) * trans(Q)

      int nr = M.Length; int nc = M[0].Length; // aka m, n
      if (nr "lt" nc)
        Console.WriteLine("ERROR: Works only m "gte" n");

      double[][] Q; double[][] R;
      MatDecompQR(M, out Q, out R); // Gram-Schmidt

      double[][] Rinv = MatInvUpperTri(R); // std algo
      double[][] Qinv = MatTranspose(Q);  // is inv(Q)
      double[][] result = MatProduct(Rinv, Qinv);
      return result;
    }

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

    private static void MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // QR decomposition, modified Gram-Schmidt
      // 'reduced' mode (all machine learning scenarios)
      int m = A.Length; int n = A[0].Length;
      if (m "lt" n)
        throw new Exception("m must be gte n ");

      double[][] QQ = new double[m][];  // working Q mxn
      for (int i = 0; i "lt" m; ++i)
        QQ[i] = new double[n];

      double[][] RR = new double[n][];  // working R nxn
      for (int i = 0; i "lt" n; ++i)
        RR[i] = new double[n];

      for (int k = 0; k "lt" n; ++k) // main loop each col
      {
        double[] v = new double[m];
        for (int i = 0; i "lt" m; ++i)  // col k
          v[i] = A[i][k];
        for (int j = 0; j "lt" k; ++j)  // inner loop
        {
          double[] colj = new double[QQ.Length];
          for (int i = 0; i "lt" colj.Length; ++i)
            colj[i] = QQ[i][j];

          double vecdot = 0.0;
          for (int i = 0; i "lt" colj.Length; ++i)
            vecdot += colj[i] * v[i];
          RR[j][k] = vecdot;

          // v = v - (R[j, k] * Q[:, j])
          for (int i = 0; i "lt" v.Length; ++i)
            v[i] = v[i] - (RR[j][k] * QQ[i][j]);
        } // j

        double normv = 0.0;
        for (int i = 0; i "lt" v.Length; ++i)
          normv += v[i] * v[i];
        normv = Math.Sqrt(normv);
        RR[k][k] = normv;

        // Q[:, k] = v / R[k, k]
        for (int i = 0; i "lt" QQ.Length; ++i)
          QQ[i][k] = v[i] / (RR[k][k] + 1.0e-12);

      } // k

      Q = QQ;
      R = RR;
    }

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

    private static double[][] MatInvUpperTri(double[][] A)
    {
      // used to invert R from QR
      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[][] 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[][] 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[][] 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 QRGramSchmidt

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

} // ns

The bottom line is: If I’m implementing linear regression from scratch and the data has categorical predictors, how the data is encoded depends on whether training will be done using stochastic gradient descent (SGD), or will be done using a matrix inverse technique, typically relaxed Moore-Penrose pseudo-inverse, or left pseudo-inverse via normal equations.

For SGD training, categorical data can be one-hot encoded. For inverse-based training techniques, data should be drop-first encoded (also called dummy encoding).

As is often the case, explaining what a complex problem is, is far more difficult than describing the answer. So bear with me.

Suppose you have some raw data that looks like:

F,short,24,arkansas,29500,liberal
M,tall,39,delaware,51200,moderate
F,short,63,colorado,75800,conservative
M,medium,36,illinois,44500,moderate
F,short,27,colorado,28600,liberal
. . .

Each line represents a person. The fields are sex, height, age, State (only Arkansas, Colorado, Delaware, Illinois), income, and political leaning. The goal is to predict income from the other variables, and you need to use linear regression.

The categorical variables, sex, height, politics, need to be encoded. The numeric predictor and the target y value should be normalized. The most common way to encode the categorical height, State, and politics variables is to use one-hot encoding. The most common way to encode the Boolean sex variable is to use 0-1 encoding. If you use divide-by-constant normalization (100 for age, 100,000 for income), the resulting normalized and 0-1, one-hot encoded data looks like:

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

OK, so far so good. But as it turns out, this normalization and encoding scheme works for SGD training but often fails for any training technique that involves a matrix inverse.

By “fails”, I mean the training will either blow up entirely, or succeed but give model weights and bias values with extremely large magnitudes, not close to the the values obtained by other training techniques.

SGD training is iterative. You specify a learning rate and a number of training epochs, and the model weights and the bias get closer and closer to optimal values. The downside here is the learning rate and max epochs values must be determined by trial and error. But SGD works with the one-hot encoding scheme.

The pseudo-inverse via normal equations technique, expressed in an equation is:

w = inv(Xt * X) * Xt * y

The w is a vector of the bias (in cell [0]) and weights (the remaining cells) that you want. X is a design matrix of predictor values, with a leading column of 1.0 values added to account for the model bias. Xt is the transpose of the design matrix. The inv() is a standard matrix inverse, but because Xt * X is mathematically positive-definite if you add a small regularization constant such as 1.0e-8 to the diagonal, you can use the Cholesky decomposition matrix inverse, which is much less complicated than alternatives. The first two * are matrix multiplication. The last * is matrix-to-vector multiplication.

It turns out that a one-hot encoded matrix of predictor values, X, is almost always ill-conditioned meaning you can’t reliably compute the inverse of X. (This last sentence is an enormous topic that would take pages to explain thoroughly). But pseudo-inverse via normal equations technique sometimes works with one-hot encoded predictors mostly Xt * X is square, and when you add a regularization constant to the diagonal, you also condition the matrix so that Cholesky inv() works.

The standard pseudo-inverse technique does not work. Expressed in an equation:

w = pinv(X) * y

The pinv() is relaxed Moore-Penrose pseudo-inverse. There are over a dozen variations of MP pinv(), all of them blisteringly complicated. Examples include SVD via Golub-Reinsch, SVD via Jacobi, SVD via QR Householder, QR via modified Gram-Schmidt, QR via Givens, QR via Householder. The standard SVD implementation in LAPACK is over 10,000 lines of Fortran code.

If you implement a MP pinv() function from-scratch, like I often do and use anything other than the insanely complicated SVD, the problem is that the design matrix X is not square and there is no way to condition the matrix, and pinv() will usually fail for any algorithm other than a most robust SVD algorithm which requires a minimum of 1,000 lines of very complicated code.

To recap to this point, for linear regression implemented from scratch, if you one-hot encode categorical predictors, SGD training works, pseudo-inverse via normal equations training sometimes works, but standard Moore-Penrose pseudo-inverse usually fails.

The most common way to deal with an ill-conditioned matrix that had one-hot encoded values, is to simply drop the first one of each encoding. For example, for variable Height, basic one-hot encoding is short = 1 0 0, medium = 0 1 0, tall = 0 0 1. If you drop the first categorical value, short, you get medium = 1 0, tall = 0 1, and short is the reference value and is assigned short = 0 0.

Dropping the first of each one-hot encoded categorical value removes the collinearity and allows matrix inverse and pseudo-inverse to work. So, if the raw data is:

F,short,24,arkansas,29500,liberal
M,tall,39,delaware,51200,moderate
F,short,63,colorado,75800,conservative
M,medium,36,illinois,44500,moderate
F,short,27,colorado,28600,liberal
. . .

Then if encoded using one-hot for SGD training:

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

But if encoded using drop-first for relaxed MP pseudo-inverse or pseudo-inverse via normal equations:

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

Interesting stuff.

I know a lot about linear regression. I know next to nothing about Art. But I can get a gut instinct about images that I think are nice in some subjective way. I like black and white shadow photography when it’s done well. But based on many of the images I see on the Internet, it’s more difficult than I expected to get an attractive image. Here are three images I like that have a kind of linear theme.

I couldn’t sleep one recent evening. So I figured I’d implement a matrix pseudo-inverse function using C#, via a recent implementation of singular value decomposition function using the Householder + QR technique. Bottom line: My implementation didn’t work as well as my older SVD using the Jacobi method.

In ordinary arithmetic, the inverse of 4 is 0.25 because 4 * 0.25 = 1. In matrix arithmetic, the inverse of a square n-by-n matrix A is an n-by-n matrix Ainv where A * Ainv = I. (Also Ainv * A = I). The * is matrix multiplication, the I is the n-by-n identity matrix which has 1s on the upper-left to lower-right diagonal, 0s elsewhere.

But for matrices that aren’t square, such as a 3-by-5 or a 4-by-2 matrix, there’s no inverse. But it’s sometimes useful to define a pseudo-inverse. There are several different types of pseudo-inverse, but the most common is called the Moore-Penrose pseudo-inverse.

For a matrix A, the MP pseudo-inverse is a matrix Apinv where A * Apinv * A = A. Computing a pseudo-inverse is tricky. One technique I sometimes prefer is the singular value decomposition (SVD) method. SVD breaks a matrix A down into three components: a matrix U, a matrix Vh, and a vector s. If you compute U * S * Vh you get the original matrix A. Here S is a square matrix with the elements of vector s on the diagonal, 0s elsewhere. The values in the s vector are called the singular values of the matrix.

Computing an SVD decomposition is very difficult. Based on my experience (but no hard data), the two most common techniques are the Jacobi method and the QR decomposition method. My recent SVD implementation uses the QR technique. (There are three common versions of QR — Householder, modified Gram-Schmidt, Givens. I used Householder which is theoretically the most stable of the three).

If you can compute the U, Vh, and S from a source matrix A, then the pseudo-inverse of A is V * Sinv * Uh. Here V is the inverse of Vh, which fortunately is just the transpose (rows and columns switched) of Vh. The Uh is the transpose of U, which again, fortunately is just the transpose of U. The Sinv is the S matrix but each value on the diagonal is 1/s instead of s.

To summarize, in order to find the MP pseudo-inverse of a matrix A, perform SVD decomposition to get U, Vh, and s. Compute V = transpose(Vh), Uh = transpose(U), Sinv (as described above), and then compute Apinv = V * Sinv * Uh.

This code shows how easy it is if you have an SVD decomposition function:

public static double[][] MatPinv(double[][] A)
{
  double[][] U; double[] s; double[][] Vh;
  MatDecompSVD(A, out U, out s, out Vh);

  double[][] Uh = MatTranspose(U);
  double[][] V = MatTranspose(Vh);
  double[][] S = MatMake(s.Length, s.Length);
  for (int i = 0; i < s.Length; ++i)
    S[i][i] = 1.0 / (s[i] + 1.0e-12); // avoid div by 0
  double[][] result = MatProduct(V, MatProduct(S, Uh));
  return result;
}

Unfortunately, computing the SVD decomposition is very, very difficult. Luckily, as I mentioned above, I had previously implemented SVD using the Householder + QR technique and so all the hard work was done.

Here’s the output of my demo program:

Begin pseudo-inverse via SVD (QR) using C#

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

Computing Moore-Penrose pseudo-inverse
Done

Pseudo-inverse:
   0.091074  -0.056097   0.165008  -0.025042  -0.014424   0.144469
   0.092124  -0.041739   0.075457  -0.137851   0.093893   0.005337
  -0.161758   0.043008  -0.142053   0.104065   0.078951  -0.041449
   0.066366  -0.025042  -0.014465   0.075336  -0.037227  -0.045420
   0.180763   0.037573   0.062447  -0.054091  -0.047030   0.010359

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

Source (wide) matrix:
   3.0   1.0   1.0   0.0   5.0
  -1.0   3.0   1.0  -2.0   4.0
   0.0   2.0   2.0   1.0  -3.0

Computing Moore-Penrose pseudo-inverse
Done

Pseudo-inverse:
   0.190177  -0.148152   0.066835
   0.001620   0.125468   0.153924
   0.064810   0.018734   0.156962
   0.084962  -0.108253   0.071392
   0.072608   0.060051  -0.102278

End demo

I validated the results using the NumPy library np.linalg.pinv() function. Note: There are subtle differences between computing a true mathematically valid Moore-Penrose pseudo-inverse and computing a "sort-of Moore-Penrose pseudo-inverse that's good enough for machine learning scenarios". I ran some tests using this pseudo-inverse code and they mostly passed, but there are still bugs, so this version is not good enough for anything other than experimentation.

An interesting way to spend time during a sleepless evening.

"The Pike" was an amusement zone in Long Beach, California, from 1902 to 1979. I went there several times as a young man. There was an interesting game called "Looff's Lite-a-Line" at the Pike. There were 15 seats for 15 players. The goal was to shoot a pinball multiple times to get five-in-a-row on a bingo-like score card before any other player did so. I think experiences like playing this game influenced my love of mathematics, which led to my love of computer science, which led to this blog site.

Demo program. Very long, very complex. Replace "lt" (less than), "gt", "lte", "gte" with Boolean operator symbols (my blog editor often chokes on symbols).

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

// no significant .NET dependencies

namespace MatrixPseudoInverseSingularValueDecompQR
{
  internal class PinvSVDProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin pseudo-inverse" +
        " via SVD (QR) using C# ");
      double[][] A = new double[6][];
      A[0] = new double[] { 1, 2, 3, 4, 5 };
      A[1] = new double[] { 0, -3, 5, -7, 9 };
      A[2] = new double[] { 2, 0, -2, 0, -2 };
      A[3] = new double[] { 4, -1, 5, 6, 1 };
      A[4] = new double[] { 3, 6, 8, 2, 2 };
      A[5] = new double[] { 5, -2, 4, -4, 3 };

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

      Console.WriteLine("\nComputing Moore-Penrose" +
        " pseudo-inverse ");
      double[][] Apinv = SVD.MatPinv(A);
      Console.WriteLine("Done ");

      Console.WriteLine("\nPseudo-inverse: ");
      MatShow(Apinv, 6, 11);

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

      double[][] B = new double[3][];
      B[0] = new double[] { 3, 1, 1, 0, 5 };
      B[1] = new double[] { -1, 3, 1, -2, 4 };
      B[2] = new double[] { 0, 2, 2, 1, -3 };
      
      Console.WriteLine("\nSource (wide) matrix: ");
      MatShow(B, 1, 6);

      Console.WriteLine("\nComputing Moore-Penrose" +
        " pseudo-inverse ");
      double[][] Bpinv = SVD.MatPinv(B);
      Console.WriteLine("Done ");

      Console.WriteLine("\nPseudo-inverse: ");
      MatShow(Bpinv, 6, 11);

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

    // helpers for Main()

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

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

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

    static double[][] MatLoad(string fn, int[] usecols,
      char sep, string comment)
    {
      // ex: double[][] A =
      // MatLoad("data.txt", new int[] { 0,1,2,3 }, ',', "#");
      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();
    }

  } // class Program

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

  public class SVD
  {
    public static double[][] MatPinv(double[][] A)
    {
      double[][] U; double[] s; double[][] Vh;
      MatDecompSVD(A, out U, out s, out Vh);

      double[][] Uh = MatTranspose(U);
      double[][] V = MatTranspose(Vh);
      double[][] S = MatMake(s.Length, s.Length);
      for (int i = 0; i "lt" s.Length; ++i)
        S[i][i] = 1.0 / (s[i] + 1.0e-12); // avoid div by 0
      double[][] result = MatProduct(V, MatProduct(S, Uh));
      return result;
    }

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

    public static void MatDecompSVD(double[][] A,
      out double[][] U, out double[] s, out double[][] Vt,
      double tol = 1.0e-12)
    {
      int m = A.Length; int n = A[0].Length;
      if (m "gte" n)
      {
        double[][] Q; double[][] R;
        MatDecompQR(A, out Q, out R);

        double[][] RtR = MatProduct(MatTranspose(R), R);

        double[] evals; double[][] V;
        MatDecompEigen(RtR, out evals, out V);

        int[] sortIdxs = VecReversed(ArgSort(evals));
        evals = VecSortedUsingIdxs(evals, sortIdxs);
        V = MatSortedUsingIdxs(V, sortIdxs);

        double[] ss = VecSqrt(VecClip(evals, 0.0));
        int[] idxs = VecIndicesWhereNonNegative(ss, tol);
        double[] ssNonZero = VecRemoveZeros(ss, idxs);
        double[][] VNonZero = MatRemoveCols(V, idxs);

        double[][] tmp1 = MatProduct(R, VNonZero);
        double[][] tmp2 = MatDivVector(tmp1, ssNonZero);
        double[][] UU = MatProduct(Q, tmp2);

        U = UU;
        s = ssNonZero;
        Vt = MatTranspose(VNonZero);
      } // m "gte" n
      else
      {
        // use the transpose trick
        double[][] dummyU;
        double[] dummyS;
        double[][] dummyVt;
        MatDecompSVD(MatTranspose(A),
          out dummyU, out dummyS, out dummyVt);
        U = MatTranspose(dummyVt);
        s = dummyS;
        Vt = MatTranspose(dummyU);
      } // m "lt" n

    } // MatDecomp()

    // ------------------------------------------------------
    // primary helpers: MatDecompQR() and MatDecompEigen()
    // ------------------------------------------------------

    private static void MatDecompQR(double[][] A,
      out double[][] Q, out double[][] R)
    {
      // reduced QR decomp using Householder reflections
      int m = A.Length; int n = A[0].Length;
      double[][] QQ = MatIdentity(m);  // working Q
      double[][] RR = MatCopy(A);  // working R

      for (int k = 0; k "lt" n; ++k)
      {
        double[] x = MatColToEnd(RR, k); // RR[k:, k]
        double normx = VecNorm(x);
        if (Math.Abs(normx) "lt" 1.0e-12) continue;

        double[] v = VecCopy(x);
        double sign;
        if (x[0] "gte" 0.0) sign = 1.0; else sign = -1.0;
        v[0] += sign * normx;
        double normv = VecNorm(v);
        for (int i = 0; i "lt" v.Length; ++i)
          v[i] /= normv;

        // apply reflection to R
        double[][] tmp1 = 
          MatRowColToEnd(RR, k); // RR[k:, k:]
        double[] tmp2 = VecMatProduct(v, tmp1);
        double[][] tmp3 = VecOuter(v, tmp2);
        double[][] B = MatScalarMult(2.0, tmp3);
        MatSubtractCorner(RR, k, B);  // R[k:, k:] -= B

        // accumulate Q
        double[][] tmp4 = 
          MatExtractAllRowsColToEnd(QQ, k); // QQ[:, k:]
        double[] tmp5 = MatVecProduct(tmp4, v);
        double[][] tmp6 = VecOuter(tmp5, v);
        double[][] C = MatScalarMult(2.0, tmp6);
        MatSubtractRows(QQ, k, C); // QQ[:, k:] -= C
      } // for k

      // return parts of QQ and RR
      Q = MatExtractFirstCols(QQ, n); // Q[:, :n]
      R = MatExtractFirstRows(RR, n); // R[:n, :]
    } // MatDecompQR()

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

    private static void MatDecompEigen(double[][] A,
      out double[] evals, out double[][] Evecs,
      double tol = 1.0e-12, int maxIter = 5000)
    {
      // eigen-decomp of symmetric matrix using QR
      int m = A.Length; int n = A[0].Length;
      double[][] AA = MatCopy(A);  // don't destroy A
      double[][] V = MatIdentity(n);

      int iter = 0;
      while (iter "lt" maxIter)
      {
        double[][] Q;
        double[][] R;
        MatDecompQR(AA, out Q, out R);
        double[][] Anew = MatProduct(R, Q);
        V = MatProduct(V, Q);

        // # check convergence (off-diagonal norm)
        double[] tmp1 = MatDiag(Anew);
        double[][] tmp2 = MatCreateDiag(tmp1);
        double[][] tmp3 = MatSubtract(Anew, tmp2);
        double norm = MatNorm(tmp3);
        if (norm "lt" tol)
        {
          AA = Anew; break;
        }
        AA = Anew;
        ++iter;
      } // iter

      if (iter == maxIter)
        Console.WriteLine("Warn: exceeded maxIter ");

      evals = MatDiag(AA);
      Evecs = V;
    } // MatDecompEigen()

    // ------------------------------------------------------
    // 32 secondary helpers
    // ------------------------------------------------------

    private static int[] ArgSort(double[] vec)
    {
      int n = vec.Length;
      int[] idxs = new int[n];
      for (int i = 0; i "lt" n; ++i)
        idxs[i] = i;

      double[] dup = new double[n];
      for (int i = 0; i "lt" n; ++i)
        dup[i] = vec[i]; // don't modify vec

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

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

    private static double[][] MatMake(int nr, int nc)
    {
      double[][] result = new double[nr][];
      for (int i = 0; i "lt" nr; ++i)
        result[i] = new double[nc];
      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[][] MatCopy(double[][] A)
    {
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[i][j];
      return result;
    }

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

    private static double[] MatDiag(double[][] A)
    {
      // return diag elements of A as a vector
      int m = A.Length; int n = A[0].Length; // m == n
      double[] result = new double[m];
      for (int i = 0; i "lt" m; ++i)
        result[i] = A[i][i];
      return result;
    }

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

    private static double[][] MatCreateDiag(double[] v)
    {
      // create a matrix using v as diagonal elements
      int n = v.Length;
      double[][] result = MatMake(n, n);
      for (int i = 0;i "lt" n; ++i)
        result[i][i] = v[i];
      return result;
    }

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

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

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

    private static double[][] MatSubtract(double[][] A,
      double[][] B)
    {
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[i][j] - B[i][j];
      return result;
    }

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

    private static double[] MatColToEnd(double[][] A,
      int k)
    {
      // last part of col k, from [k,k] to end
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[m - k];
      for (int i = 0; i "lt" m - k; ++i)
        result[i] = A[i + k][k];
      return result;
    }

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

    private static double[][] MatRowColToEnd(double[][] A,
      int k)
    {
      // block of A, row k to end, col k to end
      // A[k:, k:]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m-k, n-k);
      for (int i = 0; i "lt" m - k; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          result[i][j] = A[i + k][j + k];
      return result;
    }

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

    private static double[] MatVecProduct(double[][] A,
      double[] v)
    {
      // A * v
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[m];
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i] += A[i][j] * v[j];
      return result;
    }

    private static double[] VecMatProduct(double[] v,
      double[][] A)
    {
      // v * A
      int m = A.Length; int n = A[0].Length;
      double[] result = new double[n];
      for (int j = 0; j "lt" n; ++j)
        for (int i = 0; i "lt" m; ++i)
          result[j] += v[i] * A[i][j];
      return result;
    }

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

    private static double[][] MatScalarMult(double x,
      double[][] A)
    {
      // x * A element-wise
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = x * A[i][j];
      return result;
    }

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

    private static void MatSubtractCorner(double[][] A,
      int k, double[][] C)
    {
      // subtract elements in C from lower right A
      // A[k:, k:] -= C
      int m = A.Length; int n = A[0].Length;
      for (int i = 0; i "lt" m - k; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          A[i + k][j + k] -= C[i][j];

      return;  // modify A in-place 
    }

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

    private static void MatSubtractRows(double[][] A,
      int k, double[][] C)
    {
      // A[:, k:] -= C
      int m = A.Length; int n = A[0].Length;
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n-k; ++j)
          A[i][j+k] -= C[i][j];

      return;  // modify A in-place
    }

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

    private static double[][] 
      MatExtractAllRowsColToEnd(double[][] A, int k)
    {
      // A[:, k:]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n - k);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n - k; ++j)
          result[i][j] = A[i][j + k];
      return result;
    }

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

    private static double[][] 
      MatExtractFirstCols(double[][] A, int k)
    {
      // A[:, :n]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, k);
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" k; ++j)
          result[i][j] = A[i][j];
      return result;
    }

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

    private static double[][] 
      MatExtractFirstRows(double[][] A, int k)
    {
      // A[:n, :]
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(k, n);
      for (int i = 0; i "lt" k; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = A[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;
    }

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

    private static double[][] MatRemoveCols(double[][] A,
      int[] idxs)
    {
      // remove the columns specified in idxs
      int m = A.Length; int n = A[0].Length;
      int countNonzeos = 0;
      for (int i = 0; i "lt" idxs.Length; ++i)
        if (idxs[i] == 1) ++countNonzeos;

      double[][] result = MatMake(m, countNonzeos);
      int k = 0; // destination col
      for (int j = 0; j "lt" n; ++j)
      {
        if (idxs[j] == 0) continue; // skip this col
        for (int i = 0; i "lt" m; ++i)
        {
          result[i][k] = A[i][j];
        }
        ++k;
      }
      return result;
    }

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

    private static double[][] MatDivVector(double[][] A,
      double[] v)
    {
      // v must be all non-zero
      // divide each row of A by v, element-wise
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
        for (int j= 0; j "lt" n; ++j)
          result[i][j] = A[i][j] / v[j];
      return result;
    }

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

    private static double MatNorm(double[][] A)
    {
      // treat A as one big vector
      int m = A.Length; int n = A[0].Length;
      double sum = 0.0;
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          sum += A[i][j] * A[i][j];
      return Math.Sqrt(sum); 
    }

    private static double VecNorm(double[] v)
    {
      // standard vector norm
      double sum = 0.0;
      for (int i = 0; i "lt" v.Length; ++i)
        sum += v[i] * v[i];
      return Math.Sqrt(sum);
    }

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

    private static double[][] 
      MatSortedUsingIdxs(double[][] A, int[] idxs)
    {
      // arrange columns of A using order in idxs
      int m = A.Length; int n = A[0].Length;
      double[][] result = MatMake(m, n);
      for (int i = 0; i "lt" m; ++i)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          int srcCol = idxs[j];
          result[i][j] = A[i][srcCol];
        }

      }
      return result;
    }

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

    private static double[] 
      VecSortedUsingIdxs(double[] v, int[] idxs)
    {
      // arrange values in v using order specified in idxs
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = v[idxs[i]];
      return result;
    }

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

    private static double[] VecCopy(double[] v)
    {
      double[] result = new double[v.Length];
      for (int i = 0; i "lt" v.Length; ++i)
        result[i] = v[i];
      return result;
    }

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

    private static double[][] VecOuter(double[] v1,
      double[] v2)
    {
      // vector outer product
      double[][] result = MatMake(v1.Length, v2.Length);
      for (int i = 0; i "lt" v1.Length; ++i)
        for (int j = 0; j "lt" v2.Length; ++j)
          result[i][j] = v1[i] * v2[j];
      return result;
    }

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

    private static int[] VecReversed(int[] v)
    {
      // return vector of v elements in reversed order
      int n = v.Length;
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = v[n-1-i];
      return result;
    }

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

    private static double[] VecClip(double[] v, double minVal)
    {
      // used to make sure no negative values
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
      {
        if (v[i] "lt" minVal)
          result[i] = minVal;
        else
          result[i] = v[i];
      }
      return result;
    }

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

    private static double[] VecSqrt(double[] v)
    {
      // sqrt each element in v
      // assumes v has no negative values
      int n = v.Length;
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = Math.Sqrt(v[i]);
      return result;
    }

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

    private static int[] 
      VecIndicesWhereNonNegative(double[] v, double tol)
    {
      // indices of non-negative values (within tolerance)
      int n = v.Length;
      int[] result = new int[n];
      for (int i = 0; i "lt" n; ++i)
      {
        if (Math.Abs(v[i]) "gt" tol)
          result[i] = 1;
      }
      return result;
    }

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

    private static double[] VecRemoveZeros(double[] v,
      int[] idxs)
    {
      // if idx = 1, val is non-zero, if 0 val is near zero
      int n = v.Length;
      List"lt"double"gt" lst = new List"lt"double"gt"();
      for (int i = 0; i "lt" n; ++i)
      {
        if (idxs[i] == 1)
          lst.Add(v[i]);
      }
      return lst.ToArray();
    }

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

  } // class SVD

} // ns

In machine learning, there are two common scenarios where you want to compute the inverse of a square symmetric positive definite matrix. (Note: a positive definite matrix is always square and symmetric so “square symmetric positive definite” is a bit redundant terminology).

The first scenario is computing the pseudo-inverse using normal equations in the equation pinv(X) = inv(Xt * X) * Xt where X is a design matrix of predictors (a design matrix has a leading column of 1.0 values added to account for a bias term), and Xt is the transpose of X. If you add a small conditioning constant (like 1.0e-8) to the diagonal of (Xt * X), it is square symmetric positive definite. This scenario arises when computing the weights and bias of a linear model, typically in linear regression or quadratic regression.

The second scenario for the inverse of positive definite matrix is when computing the inverse of a kernel matrix. A kernel matrix K is a square matrix where the entries at [i,j] and [j,i] are a kernel function (usually the radial basis function, RBF) applied to training data items [i] and [j]. If you add a small conditioning constant (like 1.0e-8) to the diagonal of K, it is square symmetric positive definite. This scenario arises when finding the model weights for kernel ridge regression or the closely related Gaussian process regression.

If you have a square symmetric positive definite matrix, you can compute its inverse using Cholesky decomposition. If matrix A is square symmetric positive definite, the Cholesky decomposition gives you a square matrix L, which is lower triangular. If you compute L * Lt, where Lt is the transpose of L, you get A.

If you use Cholesky decomposition and get L and Lt, then the inverse of the source matrix A is inv(Lt) * inv(L) because:

if A = L * Lt then

inv(A) = inv(L * Lt)
       = inv(Lt) * inv(L)

The inverse of Lt is the inverse of an upper triangular matrix, which has a special, simple algorithm. Likewise, the inverse of L is the inverse of a lower triangular matrix, which has a special algorithm.

I put together a demo. The key function is:

def my_cholesky_inv(A):
  # inv(A) = inv(Lt) * inv(L)
  L = np.linalg.cholesky(A)
  Lt = np.transpose(L)
  Linv = my_inv_lower_tri(L)
  Ltinv = my_inv_upper_tri(Lt)
  result = np.matmul(Ltinv, Linv)
  return result

To validate my function, I ran the source matrix through the general purpose np.linalg.inv() function and also by using the scipy techniques with cho_solve() and cho_factor(). The demo output:

Begin inverse using Cholesky decomp

Source matrix A =
[[0.1400 0.3100 0.1700 0.2800]
 [0.3100 0.7700 0.5200 0.6600]
 [0.1700 0.5200 0.8600 0.6800]
 [0.2800 0.6600 0.6800 0.7600]]

Inverse of A using my_cholesky_inv():
[[ 790124.6258 -197532.2265  197529.0176 -296293.5019]
 [-197532.2265   49388.9715  -49380.9388   74067.4559]
 [ 197529.0176  -49380.9388   49386.1383  -74077.9017]
 [-296293.5019   74067.4559  -74077.9017  111120.4232]]

Inverse of A using np.linalg.inv():
[[ 790124.6258 -197532.2265  197529.0176 -296293.5019]
 [-197532.2265   49388.9715  -49380.9388   74067.4559]
 [ 197529.0176  -49380.9388   49386.1383  -74077.9017]
 [-296293.5019   74067.4559  -74077.9017  111120.4232]]

Inverse of A using sp.linalg.cho_factor():
[[ 790124.6258 -197532.2265  197529.0176 -296293.5019]
 [-197532.2265   49388.9715  -49380.9388   74067.4559]
 [ 197529.0176  -49380.9388   49386.1383  -74077.9017]
 [-296293.5019   74067.4559  -74077.9017  111120.4232]]

When I examined the results with precision=8 (not shown), I was very surprised at the relatively poor precision of np.linalg.inv() function on the square symmetric positive definite demo matrix — only to about 6 decimals. The moral is that if you have a square symmetric positive definite matrix, to compute the inverse you should use Cholesky decomposition inverse instead of a general purpose matrix inverse function.

One of the nice characteristics of the NumPy and SciPy libraries is that they have a consistent interface that allows various functions to be combined and chained together.

When my son was a young boy, he loved to play with his Thomas the Tank Engine set. A lot of adult men take those sets quite seriously. A fellow from the Pacific Northwest, Tom Stephenson, created a Web site (wtrak.org) where he describes a set of track modules that can be connected.

Demo program:

# cholesky_inverse_demo.py

# inverse of a (square symmetric) positive definite matrix
# using Cholesky decomposition

import numpy as np

def my_cholesky_inv(A):
  # inv(A) = inv(Lt) * inv(L)
  L = np.linalg.cholesky(A)
  Lt = np.transpose(L)
  Linv = my_inv_lower_tri(L)
  Ltinv = my_inv_upper_tri(Lt)
  result = np.matmul(Ltinv, Linv)
  return result
  
def my_inv_lower_tri(A):
  # inverse of a lower triangular matrix
  n = len(A)
  result = np.eye(n)  # Identity matrix
  for k in range(n):
    for j in range(n):
      for i in range(k):
        result[k,j] -= result[i,j] * A[k,i]
      result[k,j] /= A[k,k]
  return result

def my_inv_upper_tri(A):
  # inverse of upper triangular
  n = len(A)
  result = np.eye(n)  # Identity matrix
  for k in range(n):
    for j in range(n):
      for i in range(k):
        result[j,k] -= result[j,i] * A[i,k]
      result[j,k] /= A[k,k]
  return result

print("\nBegin inverse using Cholesky decomp ")
np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed')

# make a square symmetric positive definite matrix
B = np.array([
  [0.1, 0.2, 0.3],
  [0.4, 0.6, 0.5],
  [0.9, 0.1, 0.2],
  [0.6, 0.2, 0.6]])

A = np.matmul(B, np.transpose(B))  # A is sq symm

for i in range(len(A)):  # make it PD
  A[i,i] += 1.0e-6
print("\nSource matrix A = ")
print(A)

Ainv = my_cholesky_inv(A)
print("\nInverse of A using my_cholesky_inv(): ")
print(Ainv)

Ainv = np.linalg.inv(A)
print("\nInverse of A using np.linalg.inv(): ")
print(Ainv)

# scipy has helpers
from scipy.linalg import cho_solve, cho_factor
L = cho_factor(A, lower=True)
I = np.eye(len(A))
Ainv = cho_solve(L, I)
print("\nInverse of A using sp.linalg.cho_factor(): ")
print(Ainv)

It is possible to implement almost any kind of machine learning algorithm using C# or the Python NumPy numeric library. NumPy has hundreds of built-in functions to do almost anything. But raw C# doesn’t have any and so it’s necessary to implement even the most basic functions.

I was sitting at home one evening and speculated aboout the absolutely essential matrix and functions needed for C# machine learning. I came up with 17 functions: MatMake(), MatLoad(), MatShow(), VecShow(), MatCopy(), MatIdentity(), MatTranspose(), MatProduct(), MatSum(), MatVecProduct(), VecMatProduct(), VecToMat(), MatToVec(), VecNorm(), VecDot(), VecEucDist(), Shuffle().

Of course there are dozens more, but for a beginner I think these are definitely the required ones. And note that I’m not considering sophisticated matrix functions such as matrix inverse, matrix pseudo-inverse, and matrix decomposition.

My matrix functions use a simple array-of-arrays design, as opposed to some sort of program-defined Matrix class definition, which is overkill in my opinion.

Output of a super short demo:

Begin matrix and vector basics

A =
  1.0  2.0  3.0  4.0
  5.0  6.0  7.0  8.0
  9.0  0.0  1.0  2.0

At =
  1.0  5.0  9.0
  2.0  6.0  0.0
  3.0  7.0  1.0
  4.0  8.0  2.0

At * A =
 107.0  32.0  47.0  62.0
  32.0  40.0  48.0  56.0
  47.0  48.0  59.0  70.0
  62.0  56.0  70.0  84.0

v =
   2.0   4.0   6.0   8.0

A * v =
   60.0  140.0   40.0

End demo

The purpose of most of these 17 basic functions should be clear from their names. MatLoad() loads a matrix from values in a text file. The Shuffle() function scrambles the order of indices in an array of integers — this is used to randomize the order in which rows of training data are processed for SGD training.

I enjoy trivia contests. There are a lot of areas that are essential knowledge, such as geography, American history, and classical art. James Tissot (1836-1902) makes occassional appearances in trivia answers. Tissot did a wide range of paintings, illustrations, and caricatures. Here are two of my favorite Tissot paintings. Left: “An Interesting Story” (1872). Right: “The Gallery of HMS Calcutta” (1876).

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 MatrixVectorBasics
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin matrix and vector basics ");

      //double[][] A = MatLoad("data.txt",
      //  new int[] { 0, 1, 2, 3 }, ',', "#");

      double[][] A = new double[3][];
      A[0] = new double[] { 1, 2, 3, 4 };
      A[1] = new double[] { 5, 6, 7, 8 };
      A[2] = new double[] { 9, 0, 1, 2 };

      Console.WriteLine("\nA = ");
      MatShow(A, 1, 5);

      double[][] At = MatTranspose(A);
      Console.WriteLine("\nAt = ");
      MatShow(At, 1, 5);

      double[][] AtA = MatProduct(At, A);
      Console.WriteLine("\nAt * A = ");
      MatShow(AtA, 1, 6);

      double[] v = new double[] { 2.0, 4.0, 6.0, 8.0 };
      Console.WriteLine("\nv = ");
      VecShow(v, 1, 6);

      double[] Av = MatVecProd(A, v);
      Console.WriteLine("\nA * v = ");
      VecShow(Av, 1, 7);

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

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

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

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

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

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

    static void MatShow(double[][] A, int dec, int wid)
    {
      int nRows = A.Length; int nCols = A[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 = A[i][j];
          if (Math.Abs(v) "lt" small) v = 0.0;
          Console.Write(v.ToString("F" + dec).
            PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

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

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

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

    static double[][] MatCopy(double[][] A)
    {
      int nRows = A.Length; int nCols = A[0].Length;
      double[][] result = MatMake(nRows, nCols);
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = A[i][j];
      return result;
    }

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

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

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

    static double[][] MatTranspose(double[][] A)
    {
      int nRows = A.Length; int nCols = A[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] = A[i][j];
      return result;
    }

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

    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) // could use bRows
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

    static double[][] MatSum(double[][] A, double[][] B)
    {
      int nRows = A.Length; int nCols = A[0].Length;
      double[][] result = MatMake(nRows, nCols);
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = A[i][j] + B[i][j];
      return result;

    }
    // --------------------------------------------------------

    static double[] MatVecProd(double[][] A, double[] v)
    {
      // return a regular vector
      int nRows = A.Length; int nCols = A[0].Length;
      int n = v.Length;
      if (nCols != n)
        throw new Exception("non-comform in MatVecProd");

      double[] result = new double[nRows];
      for (int i = 0; i "lt" nRows; ++i)
        for (int k = 0; k "lt" nCols; ++k)
          result[i] += A[i][k] * v[k];

      return result;
    }

    static double[] VecMatProd(double[] v, double[][] A)
    {
      // return a regular vector
      int nRows = A.Length; int nCols = A[0].Length;
      int n = v.Length;

      if (n != nRows)
        throw new Exception("non-comform in VecMatProd");

      double[] result = new double[nCols];
      for (int j = 0; j "lt" nCols; ++j)
        for (int k = 0; k "lt" n; ++k)
          result[j] += v[k] * A[k][j];

      return result;
    }

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

    static double[][] VecToMat(double[] vec,
      int nRows, int nCols)
    {
      double[][] result = MatMake(nRows, nCols);
      int k = 0;
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = vec[k++];
      return result;
    }

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

    static double[] MatToVec(double[][] A)
    {
      int nRows = A.Length; int nCols = A[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++] = A[i][j];
        }
      return result;
    }

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

    static double VecNorm(double[] vec)
    {
      int n = vec.Length;
      double sum = 0.0;
      for (int i = 0; i "lt" n; ++i)
        sum += vec[i] * vec[i];
      return Math.Sqrt(sum);
    }

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

    static double VecDot(double[] v1, double[] v2)
    {
      double result = 0.0;
      int n = v1.Length;
      for (int i = 0; i "lt" n; ++i)
        result += v1[i] * v2[i];
      return result;
    }

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

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

    }

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

    static void Shuffle(int[] indices, Random rnd)
    {
      // Fisher-Yates. last iteration not necessary
      for (int i = 0; i "lt" indices.Length; ++i)
      {
        int ri = rnd.Next(i, indices.Length);
        int tmp = indices[i];
        indices[i] = indices[ri];
        indices[ri] = tmp;
      }
    } // Shuffle

  } // class Program

} // ns