James D. McCaffrey

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

I wrote an article titled “Decision Tree Regression from Scratch Without Pointers or Recursion Using C#” in the February 2026 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2026/02/17/decision-tree-regression-from-scratch.aspx.

Decision tree regression is a machine learning technique that incorporates a set of if-then rules in a tree data structure to predict a single numeric value. For example, a decision tree regression model prediction might be, “If employee age is greater than 39.0 and age is less than or equal to 42.5 and years-experience is less than or equal to 8.0 and height is greater than 64.0 then bank account balance is $754.14.”

The article presents decision tree regression, implemented from scratch with C#, using list storage for tree nodes (no pointers/references), a FIFO stack to build the tree (no recursion), weighted variance minimization for the node split function, and saves associated row indices in each node. This design has maximum flexibility.

The output of the demo program is:

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

Setting maxDepth = 3
Setting minSamples = 2
Setting minLeaf = 18
Using default numSplitCols = -1
Creating and training tree
Done

Tree:
ID 0   | idx   0 | v  -0.2102 | L   1 | R   2 | y 0.3493 | leaf F
ID 1   | idx   4 | v   0.1431 | L   3 | R   4 | y 0.5345 | leaf F
ID 2   | idx   0 | v   0.3915 | L   5 | R   6 | y 0.2382 | leaf F
ID 3   | idx   0 | v  -0.6553 | L   7 | R   8 | y 0.6358 | leaf F
ID 4   | idx  -1 | v   0.0000 | L   9 | R  10 | y 0.4123 | leaf T
ID 5   | idx   4 | v  -0.2987 | L  11 | R  12 | y 0.3032 | leaf F
ID 6   | idx   2 | v   0.3777 | L  13 | R  14 | y 0.1701 | leaf F
ID 7   | idx  -1 | v   0.0000 | L  15 | R  16 | y 0.6952 | leaf T
ID 8   | idx  -1 | v   0.0000 | L  17 | R  18 | y 0.5598 | leaf T
ID 11  | idx  -1 | v   0.0000 | L  23 | R  24 | y 0.4101 | leaf T
ID 12  | idx  -1 | v   0.0000 | L  25 | R  26 | y 0.2613 | leaf T
ID 13  | idx  -1 | v   0.0000 | L  27 | R  28 | y 0.1882 | leaf T
ID 14  | idx  -1 | v   0.0000 | L  29 | R  30 | y 0.1381 | leaf T

Evaluating model

Accuracy train (within 0.10) = 0.3750
Accuracy test (within 0.10) = 0.4750

MSE train = 0.0048
MSE test = 0.0054

Predicting for trainX[0] =
  -0.1660   0.4406  -0.9998  -0.3953  -0.7065
Predicted y = 0.4101

IF
column 0  "gt "   -0.2102 AND
column 0  "lte"   0.3915 AND
column 4  "lte"  -0.2987 AND
THEN
predicted = 0.4101

One advantage of decision tree regression compared to other techniques, such as kernel ridge regression and neural network regression, is that decision trees are significantly more interpretable. This diagram shows how the tree is structured and how the prediction is made:

A single decision tree can be used by itself for regression problems. But a more common approach is to use a collection of decision trees, called an ensemble. Examples include bagging tree regression, random forest regression, adaptive boosting regression, and gradient boosting regression.

Left: “Maleficent” (2014)
Center: “A Monster Calls” (2016)
Right: “The Lord of the Rings: The Two Towers” (2002)

Bottom line: I recently implemented a function to compute the singular value decomposition (SVD) of a matrix, from scratch, using Python. That implementation used many built-in NumPy functions such as np.diag(). I refactored the implementation to remove all those NumPy calls, making a version of SVD that’s even more from-scratch. This implementation uses the Householder + QR method. It works but it’s not as stable as my older version that uses the Jacobi algorithm. Put another way, the code in this blog post is interesting but not practical.

Computing the singular value decomposition (SVD) of a matrix is one of the most difficult problems in numerical programming. Many researchers worked for many years to develop a solution. The standard version of SVD (LAPACK) is over 10,000 lines of Fortran code!

I’ve plunked away at a from-scratch implementation of SVD for many years. When I got some free time recently, I was able to dedicate about 50 hours of my time for a pretty decent (estimated accuracy to about 1.0e-10 and reasonable performance on matrices up to about 1000 rows).

The first part of my demo shows an example of SVD on a tall (more rows than columns) matrix:

Begin SVD demo

TALL MATRIX

A =
[[ 1  2  3  4  5]
 [ 0 -3  5 -7  9]
 [ 2  0 -2  0 -2]
 [ 4 -1  5  6  1]
 [ 3  6  8  2  2]
 [ 5 -2  4 -4  3]]

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]

Vt =
[[ 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]]

reconstructed A
[[ 1.000000  2.000000  3.000000  4.000000  5.000000]
 [ 0.000000 -3.000000  5.000000 -7.000000  9.000000]
 [ 2.000000  0.000000 -2.000000  0.000000 -2.000000]
 [ 4.000000 -1.000000  5.000000  6.000000  1.000000]
 [ 3.000000  6.000000  8.000000  2.000000  2.000000]
 [ 5.000000 -2.000000  4.000000 -4.000000  3.000000]]

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

np.linalg.svd results:

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]

Vt =
[[-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]]

I validated the demo by using the np.linalg.svd() function to make sure I got the same results. Note that the sign of values in the columns of U and the rows of Vt are arbitrary for SVD.

For the second part of the demo, I used a wide (more columns than rows) matrix:

WIDE MATRIX

A =
[[ 3.000000  1.000000  1.000000  0.000000  5.000000]
 [-1.000000  3.000000  1.000000 -2.000000  4.000000]
 [ 0.000000  2.000000  2.000000  1.000000 -3.000000]]

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]

Vt =
[[-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]]

reconstructed A
[[ 3.000000  1.000000  1.000000 -0.000000  5.000000]
 [-1.000000  3.000000  1.000000 -2.000000  4.000000]
 [ 0.000000  2.000000  2.000000  1.000000 -3.000000]]

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

np.linalg.svd results:

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]

Vt =
[[-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]]

End demo

The refactoring took me a long time. For example, the original call to the statement:

R[k:, k:] -= 2 * np.outer(v, v @ R[k:, k:])

uses complex slicing, matrix multiplication, the np.outer() function, and scalar multiplication. It was replaced by:

tmp1 = my_row_col_to_end(R, k)  # R[k:, k:]
tmp2 = my_vec_mat_product(v, tmp1)  # v @ R[k:, k:]
tmp3 = my_outer(v, tmp2)  # np.outer(v, v @ R[k:, k:])
B = my_scalar_mult(2.0, tmp3)
my_subtract_corner(R, k, B)  # modify R in-plce

Five statements that call five from-scratch functions.

I enjoyed the effort and learned a lot about SVD decomposition.

While I was refactoring my SVD code, I was faced with many situations where I had to decide to truncate code (to reduce number of lines) or not truncate (better clarity). Human limb truncation is usually not optional.

When I was growing up, my father was a doctor and my mother was a nurse. So my brother and sister and I grew up in a kind of medical environment. We learned early on that physical anomalies were nothing to be afraid of. Here are two amusing examples. Left: Clever math. Right: Clever tattoo.

Demo code. Caution: Almost certainly has a few minor bugs for edge cases. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor chokes on symbols).

# svd_from_scratchier.py

import numpy as np
import math

# -----------------------------------------------------------
# 31 helper functions to replace NumPy calls
# -----------------------------------------------------------

def my_norm(v):
  sum = 0.0
  for i in range(len(v)):
    sum += v[i] * v[i]
  return math.sqrt(sum)

def my_norm_matrix(A):
  m,n = A.shape
  sum = 0.0
  for i in range(m):
    for j in range(n):
      sum += A[i][j] * A[i][j]
  return math.sqrt(sum)

def my_diag(A):
  m,n = A.shape
  result = np.zeros(m)
  for i in range(m):
    result[i] = A[i][i]
  return result

def my_create_diag(v):
  # matrix from diag vector
  n = len(v)
  result = np.zeros((n,n))
  for i in range(n):
    result[i][i] = v[i]
  return result

def my_eye(n):
  result = np.zeros((n,n))
  for i in range(n):
    result[i][i] = 1.0
  return result

def my_copy_matrix(A):
  m,n = A.shape
  result = np.zeros((m,n))
  for i in range(m):
    for j in range(n):
      result[i][j] = A[i][j]
  return result 

def my_copy_vector(v):
  n = len(v)
  result = np.zeros(n)
  for i in range(n):
    result[i] = v[i]
  return result 

def my_col_to_end(A, k):
  # last part of col k, from k,k to end
  m,n = A.shape
  result = np.zeros(m-k)
  for i in range(0, m-k):
    result[i] = A[i+k][k]
  return result

def my_row_col_to_end(A, k):
  # row k to end, col k to end
  m,n = A.shape
  result = np.zeros((m-k,n-k))
  for i in range(0,m-k):
    for j in range(0,n-k):
      result[i][j] = A[i+k][j+k]
  return result

def my_vec_mat_product(v, A):
  m,n = A.shape
  result = np.zeros(n)
  for j in range(n):
    for k in range(m):
      result[j] += v[k] * A[k][j]
  return result

def my_mat_vec_product(A, v):
  m,n = A.shape
  result = np.zeros(m)
  for i in range(m):
    for k in range(n):
      result[i] +=A[i][k] * v[k]
  return result

def my_outer(v1, v2):
  result = np.zeros((v1.size, v2.size))
  for i in range(v1.size):
    for j in range(v2.size):
      result[i][j] = v1[i] * v2[j]
  return result

def my_scalar_mult(x, A):
  # x times each element A
  m,n = A.shape
  result = np.zeros((m,n))
  for i in range(m):
    for j in range(n):
      result[i][j] = x * A[i][j]
  return result 

def my_subtract_corner(A, k, C):
  # subtract elements in C from lower right A
  #  starting at [k][k]
  m,n = A.shape
  for i in range(0, m-k):
    for j in range(0, n-k):
      A[i+k][j+k] -= C[i][j]
  return  # in-place

def my_extract_all_rows_col_k_to_end(A, k):
  m,n = A.shape
  result = np.zeros((m,n-k))
  for i in range(m):
    for j in range(n-k):
      result[i][j] = A[i][j+k]
  return result

def my_subtract_all_rows_col_k_to_end(A, k, C):
  # subtract elements in C from right of A
  m,n = A.shape
  for i in range(0, m):
    for j in range(0, n-k):
      A[i][j+k] -= C[i][j] 

def my_extract_first_k_cols(A, k):
  m,n = A.shape
  result = np.zeros((m,k))
  for i in range(m):
    for j in range(k):
      result[i][j] = A[i][j]
  return result

def my_extract_first_k_rows(A, k):
  m,n = A.shape
  result = np.zeros((k,n))
  for i in range(k):
    for j in range(n):
      result[i][j] = A[i][j]
  return result

def my_mat_subtract(A, B):
  m,n = A.shape
  result = np.zeros((m,n))
  for i in range(m):
    for j in range(n):
      result[i][j] = A[i][j] - B[i][j]
  return result  
  
def my_mat_prod(A, B):
  a_rows = len(A)
  a_cols = len(A[0])
  b_rows = len(B)
  b_cols = len(B[0])
  result = np.zeros((a_rows, b_cols))
  for i in range(a_rows):
    for j in range(b_cols):
      for k in range(a_cols):
        result[i][j] += A[i][k] * B[k][j]
  return result

def my_mat_transpose(A):
  m,n = A.shape
  result = np.zeros((n,m))
  for i in range(m):
    for j in range(n):
      result[j][i] = A[i][j]
  return result

def my_argsort(v):
  n = len(v)
  augmented = []
  for i in range(n):
    tmp = [v[i],i]
    augmented.append(tmp)
  augmented.sort(key=lambda x: x[0])
  result = np.zeros(n, dtype=np.int64)
  for i in range(n):
    result[i] = augmented[i][1]
  return result 

def my_reversed(v):
  n = len(v)
  result = np.zeros(n, dtype=np.int64)
  for i in range(n):
    result[i] = v[n-1-i]
  return result 

def my_sort_using_idxs(v, idxs):
  n = len(v)
  result = np.zeros(n)
  for i in range(n):
    result[i] = v[idxs[i]]
  return result

def my_sort_cols_using_idxs(A, idxs):
  m,n = A.shape
  result = np.zeros((m,n))
  for i in range(m):
    for j in range(n):
      src_col = idxs[j]
      result[i][j] = A[i][src_col]
  return result
    
def my_clip_to_nonnegative(v):
  n = len(v)
  result = np.zeros(n)
  for i in range(n):
    if v[i] "lt" 0.0:
      result[i] = 0.0
    else:
      result[i] = v[i]
  return result

def my_sqrt(v):
  n = len(v)
  result = np.zeros(n)
  for i in range(n):
    result[i] = math.sqrt(v[i])
  return result

def my_indices_where_nonzero(v, tol):
  # 1 == not zero
  n = len(v)
  result = np.zeros(n, dtype=np.int64)
  for i in range(n):
    if v[i] "gt" tol:
      result[i] = 1  # non-zero
  return result

def my_remove_zeros(v, indices):
  # indices: 1 is non-zero, 0 is a near-zero
  n = len(v)
  lst = []
  for i in range(n):
    if indices[i] == 1:
      lst.append(v[i])
  result = np.array(lst)
  return result

def my_remove_cols(A, indices):
  m,n = A.shape
  ct_nonzeros = 0
  for i in range(len(indices)):
    if indices[i] == 1:
      ct_nonzeros += 1
  result = np.zeros((m, ct_nonzeros))

  k = 0  # destination col
  for j in range(n):
    if indices[j] == 0: continue # next src col
    for i in range(m):
      result[i][k] = A[i][j]
    k += 1 
  return result

def my_mat_div_vector(A, v):
  # A / v
  m,n = A.shape
  if n != len(v): print("FATAL in my_mat_div_vector")
  result = np.zeros((m,n))
  for i in range(m):
    for j in range(n):
      result[i][j] = A[i][j] / v[j]
  return result  

# -----------------------------------------------------------
# helpers for my_svd(): my_qr_decomp(), my_eigen()
# -----------------------------------------------------------

def my_qr_decomp(A):
  # reduced QR decomposition using Householder reflections
  # A: (m x n), m "gte" n
  # returns Q: (m x n), R: (n x n)
  # used by both my_eigen() and my_svd()

  # A = A.astype(float).copy() # assume A is float
  m, n = A.shape

  # Q = np.eye(m)
  Q = my_eye(m)
  # R = A.copy()
  R = my_copy_matrix(A)

  for k in range(n):
    # x = R[k:, k]
    x = my_col_to_end(R, k)
    # normx = np.linalg.norm(x)
    normx = my_norm(x)
    if normx == 0:
      continue

    # v = x.copy()
    v = my_copy_vector(x)
    sign = 1.0 if x[0] "gte" 0 else -1.0
    v[0] += sign * normx
    # v /= np.linalg.norm(v)
    v /= my_norm(v)

    # Apply reflector to R
    # R[k:, k:] -= 2 * np.outer(v, v @ R[k:, k:])
    tmp1 = my_row_col_to_end(R, k)  # R[k:, k:]
    tmp2 = my_vec_mat_product(v, tmp1)
    tmp3 = my_outer(v, tmp2)
    B = my_scalar_mult(2.0, tmp3)
    my_subtract_corner(R, k, B)

    # Accumulate Q
    # Q[:, k:] -= 2 * np.outer(Q[:, k:] @ v, v)
    tmp4 = my_extract_all_rows_col_k_to_end(Q, k)
    tmp5 = my_mat_vec_product(tmp4, v)
    tmp6 = my_outer(tmp5, v)
    C = my_scalar_mult(2.0, tmp6)
    my_subtract_all_rows_col_k_to_end(Q, k, C)

  # return Q[:, :n], R[:n, :]
  Q_part = my_extract_first_k_cols(Q, n)
  R_part = my_extract_first_k_rows(R, n)
  return Q_part, R_part

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

def my_eigen(A, tol=1.0e-12, max_iter=2000):
  # eigen-decomp of symmetric matrix using QR iteration
  # returns eigenvalues and eigenvectors.

  # A = A.astype(float).copy()  # assume float
  AA = my_copy_matrix(A)
  n = AA.shape[0]

  # V = np.eye(n)
  V = my_eye(n)  # eigenvectors (columns)

  for _ in range(max_iter):
    # QR decomposition (using our own QR)
    Q, R = my_qr_decomp(AA)
    # A_new = R @ Q
    A_new = my_mat_prod(R, Q)
    # V = V @ Q
    V = my_mat_prod(V, Q)

    # Check convergence (off-diagonal norm)
    # off = A_new - np.diag(np.diag(A_new))
    off = \
      my_mat_subtract(A_new, my_create_diag(my_diag(A_new)))
    # if np.linalg.norm(off) "lt" tol:
    if my_norm_matrix(off) "lt" tol:
      AA = A_new
      break

    AA = A_new

  # eigvals = np.diag(A)
  eigvals = my_diag(AA)
  return eigvals, V

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

def my_svd(A, tol=1.0e-12):
  # SVD using only NumPy operations:

  # A = np.asarray(A, float)  # assume A is float
  m, n = A.shape

  if m "gte" n:
    Q, R = my_qr_decomp(A)

    # symmetric eigenproblem
    # RtR = R.T @ R
    RtR = my_mat_prod(my_mat_transpose(R), R)
    eigvals, V = my_eigen(RtR, tol=tol)

    # sort large to small
    # idx = np.argsort(eigvals)[::-1]
    sort_idxs = my_reversed(my_argsort(eigvals))
    # eigvals = eigvals[idx]
    eigvals = my_sort_using_idxs(eigvals, sort_idxs)

    # V = V[:, idx]
    V = my_sort_cols_using_idxs(V, sort_idxs)

    # S = np.sqrt(np.clip(eigvals, 0, None))
    s = my_sqrt(my_clip_to_nonnegative(eigvals))

    # nonzero = S "gt" tol
    # U = Q @ (R @ V[:, nonzero] / S[nonzero])

    idxs = my_indices_where_nonzero(s, tol)
    s_nonzero = my_remove_zeros(s, idxs)
    V_nonzero = my_remove_cols(V, idxs)

    # U = Q @ ((R @ V_nonzero) / S_nonzero)
    # tmp1 = R @ V_nonzero
    tmp1 = my_mat_prod(R, V_nonzero)
    tmp2 = my_mat_div_vector(tmp1, s_nonzero)
    # U = Q @ tmp2
    U = my_mat_prod(Q, tmp2)

    # U = Q @ (R @ (V_nonzero / S_nonzero))

    # return U, S[nonzero], V[:, nonzero].T
    return U, s_nonzero, my_mat_transpose(V_nonzero)

  else:
    # Wide matrix: transpose trick
    # U, S, Vt = svd_qr_eig_numpy_only(A.T, tol)
    # return Vt.T, S, U.T
    U, s, Vt = my_svd(my_mat_transpose(A), tol)
    return my_mat_transpose(Vt), s, my_mat_transpose(U)
   

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

print("\nBegin SVD demo ")

print("\nTALL MATRIX ")

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

A = np.array([[3, 1, 1],
              [-1, 3, 1],
              [0, 2, 2],
              [4, 4, -2]])

A = np.array([[1, 2, 3, 4, 5],
              [0,-3, 5, -7, 9],
              [2, 0, -2, 0, -2],
              [4, -1, 5, 6, 1],
              [3, 6, 8, 2, 2],
              [5, -2, 4, -4, 3]])

print("\nA = ");print(A)

U, s, Vt = my_svd(A)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)

A_reconstructed = my_mat_prod(U, \
  my_mat_prod(my_create_diag(s), Vt))
print("\nreconstructed A"); print(A_reconstructed)

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

print("\nnp.linalg.svd results: ")
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)

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

print("\nWIDE MATRIX ")

A = np.array([[3, 1, 1, 0, 5],
              [-1, 3, 1, -2, 4],
              [0, 2, 2, 1, -3]], dtype=np.float64)

print("\nA = ");print(A)

U, s, Vt = my_svd(A)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)

A_reconstructed = my_mat_prod(U, \
  my_mat_prod(my_create_diag(s), Vt))
print("\nreconstructed A"); print(A_reconstructed)

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

print("\nnp.linalg.svd results: ")
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)

print("\nEnd demo ")

I like to write code almost every day, to entertain myself, and to keep in practice. One morning before work, I realized it had been many, many months since I implemented a function to compute the inverse of a matrix using the Gauss-Jordan elimination to reduced row echelon form.

The Wikipedia article on matrix inverse lists 14 different algorithms, and each algorithm has multiple variations, and each variation has dozens of implementation alternatives. Six examples include 1.) LUP decomposition algorithm (Crout version and others), 2.) QR decomposition algorithm (Householder version and others), 3.) SVD decomposition algorithm (Jacobi version and others), 4.) Cholesky decomposition (Banachiewicz version and others), 5.) Cayley-Hamilton algorithm (Faddeev–LeVerrier version and others), 6.) Newton Iteration. The fact that there are so many different techniques is a direct indication of how tricky it is to compute a matrix inverse.

The Gauss-Jordan technique is arguably the simplest way to compute a matrix inverse. But Gauss-Jordan is rarely used in practice because, compared to other techniques, Gauss-Jordan is susceptible to numerical instability, especially when the source matrix is nearly singular (meaning has no inverse). Gauss-Jordan inverse is mostly a teaching technique for students.

Anyway, after a few minutes of coding, I had a C# language demo up and running. The output is:

Matrix inverse using Gauss-Jordan elimination 
  to reduced row echelon form

Source matrix:
  3.0  5.0  9.0
  1.0  1.0  5.0
  2.0  4.0  8.0

Computing inverse
Done

Inverse from Gauss-Jordan:
   1.5000   0.5000  -2.0000
  -0.2500  -0.7500   0.7500
  -0.2500   0.2500   0.2500

Minv * M:
   1.0000  -0.0000  -0.0000
   0.0000   1.0000   0.0000
   0.0000   0.0000   1.0000

M * Minv:
   1.0000  -0.0000   0.0000
  -0.0000   1.0000   0.0000
  -0.0000  -0.0000   1.0000

End demo

The Gauss-Jordan technique starts by creating a matrix that is the source matrix that is augmented by the Identity matrix to the right. For the demo matrix, the augmented matrix is:

 3  5  9   1  0  0
 1  1  5   0  1  0
 2  4  8   0  0  1

Then the algorithm iterates through each row and performs three possible row operations:

* exchange two rows
* multiply a row by a constant
* multiply a row by a constant times another row

The version I implemented doesn’t exchange rows, but some versions do. The goal is to transform the augmented matrix so that the left side of the augmented matrix becomes the identity matrix. For the demo, the result is:

 1  0  0    1.50   0.50  -2.00
 0  1  0   -0.25  -0.75   0.75  
 0  0  1   -0.25   0.25   0.25

When finished, the right side of the transformed augmented matrix is the inverse of the source matrix. Clever.

OK, much fun for me. But as mentioned, the Gauss-Jordan algorithm to compute a matrix inverse is not used very often, so this was just mental exercise for me.

I loved building model airplanes when I was a young man. Two of my favorite models were World War II aircraft that featured inverted gull wing designs.

Right: The U.S. Vought F4U Corsair was a fighter bomber that entered service in mid-war. It was a rare example of a plane that was equally effective as a fighter and as an attack bomber. It was one of the fastest planes of the war with a top speed of 440 mph. The inverted gull wing design allowed a huge 13 foot 4 inch propeller.

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

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

namespace MatrixInverseGaussJordan
{
  internal class MatrixInverseGaussJordanProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nMatrix inverse " +
        "using Gauss-Jordan elimination to reduced " +
        "row echelon form ");

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


      Console.WriteLine("\nSource matrix: ");
      MatShow(M, 1, 5);

      Console.WriteLine("\nComputing inverse ");
      double[][] Minv = MatInverseGaussJordan(M);
      Console.WriteLine("Done ");

      Console.WriteLine("\nInverse from Gauss-Jordan: ");
      MatShow(Minv, 4, 9);

      double[][] Minv_M = MatMult(Minv, M);
      Console.WriteLine("\nMinv * M: ");
      MatShow(Minv_M, 4, 9);

      double[][] M_Minv = MatMult(M, Minv);
      Console.WriteLine("\nM * Minv: ");
      MatShow(M_Minv, 4, 9);

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

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

    static double[][] MatInverseGaussJordan(double[][] A)
    {
      int n = A.Length;  // assume A is square
      double[][] B = MatMake(n, 2 * n);  // make augmented
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lt" n; ++j)
        {
          B[i][j] = A[i][j];
        }
        B[i][i + n] = 1.0;
      }

      for (int i = 0; i "lt" n; ++i)
      {
        double pv = B[i][i];
        if (Math.Abs(pv) "lt" 1.0e-12)
          Console.WriteLine("FATAL: matrix is singular ");

        for (int j = 0; j "lt" 2 * n; ++j)
          B[i][j] /= pv;  // make diag element 1.0
        for (int j = 0; j "lt" n; ++j)  // elim other rows
        {
          if (i != j)
          {
            double factor = B[j][i];
            for (int k = 0; k "lt" 2 * n; ++k)
              B[j][k] -= factor * B[i][k];  // tricky
          }
        }
      } // i

      // extract the inverse from augmented
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = B[i][j + n];

      return result;
    }

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

    // helper functions:
    // MatMake(), MatShow(), MatMult(), MatLoad()

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

    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 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)
    {
      // load matrix from text file
      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[][] MatMult(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 Program

} // ns