One morning before work, I figured I’d tidy up my C# implementation of singular value decomposition (SVD). In ordinary arithmentic you can decompose the number 24 as 3 * 2 * 5. In matrix algebra, you can decompose a matrix like A = U * S * Vh.

The resulting U is a matrix, S is diagonal matrix (values on the diagonal, 0s elsewhere), and Vh is a matrix (the “h” indicates it’s the transpose of its hidden partner V).

For example, if matrix A =

  1.0000  2.0000  3.0000
  5.0000  0.0000  2.0000
  8.0000  5.0000  4.0000
  1.0000  0.0000  9.0000

then the result of SVD decomposition is:

U:
  0.2581  0.1200 -0.4902
  0.3515 -0.2318  0.8110
  0.7265 -0.5233 -0.2994
  0.5310  0.8112  0.1111

s:
 13.0781  7.1542  2.7892

Vh:
  0.6392  0.3172  0.7006
 -0.6171 -0.3322  0.7134
  0.4590 -0.8883 -0.0166

The s result is a vector. The associated S matrix is:

  13.0781  0.0000  0.0000
   0.0000  7.1542  0.0000
   0.0000  0.0000  2.7892

And if you multiply U * S * Vh you get the original matrix A back:

  1.0000  2.0000  3.0000
  5.0000  0.0000  2.0000
  8.0000  5.0000  4.0000
  1.0000  0.0000  9.0000

In machine learning scenarios, singular value decomposition doesn’t have any direct use, but SVD is used by many important ML algorithms. For example, SVD can be used to compute the Moore-Penrose matrix pseudo-inverse, which can then be used to solve for the weights of a linear regression model.

Unfortunately, SVD is one of the most complicated functions to implement in all of numerical programming. My implementation is based on one of the GNU Scientific Library versions. My implementation is the “full_matrices=False” version, as opposed to the “full_matrices=True” version. This is a very tricky topic that is outside the scope of this blog post. In general, the “full_matrices=False” version works fine for most ML scenarios (in fact, all ML scenarios in my experience).

The SVD implementation in this post is a completely standalone function, with no external dependencies and no helper functions. The implementation works only with regular numbers, not complex a+bi numbers. The calling code looks like:

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

double[][] U;
double[] s;
double[][] Vh;
MatDecompSVD(A, out U, out s, out Vh);

In matrix algebra, m usually represents the number of rows and n is the number of columns. One of the many tricky things about SVD is that implementation and results depend on whether m is greater than n (such as most ML datasets) or m is less than n (common in non-ML scenarios), or m equals n (common in basic linear algebra problems).

The demo output for an example where m is less than n is:

B:
 -1.0000  2.0000  3.0000  9.0000
  5.0000  0.0000 -2.0000  4.0000
  8.0000 -5.0000  4.0000  7.0000

U:
  0.4656  0.8850  0.0040
  0.3577 -0.1923  0.9138
  0.8095 -0.4241 -0.4061

s:
 14.6049  7.8902  4.2944

Vh:
  0.5340 -0.2134  0.2683  0.7729
 -0.6640  0.4930  0.1703  0.5358
  0.3064  0.4747 -0.8011  0.1975

U * S * Vh =
 -1.0000  2.0000  3.0000  9.0000
  5.0000  0.0000 -2.0000  4.0000
  8.0000 -5.0000  4.0000  7.0000

In addition to checking U * S * Vh, I also validated the results using the numpy linalg.svd() function, which is called like:

import numpy as np

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

(U, s, Vh) = np.linalg.svd(B, full_matrices=False)

Anyhow, good fun.

Most classic matrix algorithms were developed in the 1950s and 1960s, and many of these algorithms are elegant and beautiful. I’m not a fan of cigarette smoke, but some of the old advertising is elegant and beautiful (to my eye anyway).

The image on the right isn’t an advertisement per say. It’s the cover of Spanish Magazine “Nuevo Mundo” (“New World”) that was popular from 1894 to 1933. The art deco style cover illustrations was a showcase for many Spanish illustrators and artists.

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

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

namespace SingularDecompStandalone
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin singular value" +
        " decomposition using C# ");

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

      // 1. m (num rows) greater than n (num cols)
      double[][] A = new double[4][];
      A[0] = new double[] { 1, 2, 3 };
      A[1] = new double[] { 5, 0, 2 };
      A[2] = new double[] { 8, 5, 4 };
      A[3] = new double[] { 1, 0, 9 };

      Console.WriteLine("\nA: ");
      for (int i = 0; i "lt" A.Length; ++i)
        VecShow(A[i], 4, 8);

      double[][] U; double[] s; double[][] Vh;
      MatDecompSVD(A, out U, out s, out Vh);
      Console.WriteLine("\nU: ");
      MatShow(U, 4, 8);
      Console.WriteLine("\ns: ");
      VecShow(s, 4, 8);
      Console.WriteLine("\nVh: ");
      MatShow(Vh, 4, 8);

      double[][] S = VecToDiagMat(s);
      double[][] checkA = MatProduct(U, MatProduct(S, Vh));
      Console.WriteLine("\nU * S * Vh = ");
      MatShow(checkA, 4, 8);

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

      // 2. m less than n
      double[][] B = new double[3][];
      B[0] = new double[] { -1, 2, 3, 9 };
      B[1] = new double[] { 5, 0, -2, 4 };
      B[2] = new double[] { 8, -5, 4, 7 };

      Console.WriteLine("\nB: ");
      for (int i = 0; i "lt" B.Length; ++i)
        VecShow(B[i], 4, 8);

      // double[][] U; double[] s; double[][] Vh;
      MatDecompSVD(B, out U, out s, out Vh);
      Console.WriteLine("\nU: ");
      MatShow(U, 4, 8);
      Console.WriteLine("\ns: ");
      VecShow(s, 4, 8);
      Console.WriteLine("\nVh: ");
      MatShow(Vh, 4, 8);

      S = VecToDiagMat(s);
      double[][] checkB = MatProduct(U, MatProduct(S, Vh));
      Console.WriteLine("\nU * S * Vh = ");
      MatShow(checkB, 4, 8);

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

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();

    } // Main()

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

    static void MatDecompSVD(double[][] M,
      out double[][] U, out double[] s, out double[][] Vh)
    {
      // Jacobi algorithm
      // see github.com/ampl/gsl/blob/master/linalg/svd.c 
      // "full_matrices=False" version
      // does not sort s high to low
      double DBL_EPSILON = 1.0e-15;
      int m = M.Length;  // nRows
      int n = M[0].Length;  // nCols

      // A = MatCopy(M); 
      double[][] A = new double[m][];  // working U
      for (int i = 0; i "lt" m; ++i)
        A[i] = new double[n];
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          A[i][j] = M[i][j];

      // Q = MatIdentity(n);
      double[][] Q = new double[n][];  // working V
      for (int i = 0; i "lt" n; ++i)
        Q[i] = new double[n];
      for (int i = 0; i "lt" n; ++i)
        Q[i][i] = 1.0;
      //
      double[] t = new double[n];  // working s

      // initialize counters
      int count = 1;
      int sweep = 0;
      double tol = 10 * m * DBL_EPSILON; // heuristic

      // Always do at least 12 sweeps
      int sweepmax = Math.Max(5 * n, 12); // heuristic

      // store the column error estimates in St for use
      // during orthogonalization

      for (int j = 0; j "lt" n; ++j)
      {
        // double[] cj = MatGetColumn(A, j);
        double[] cj = new double[m];
        for (int i = 0; i "lt" m; ++i)
          cj[i] = A[i][j];

        // double sj = VecNorm(cj);
        double sj = 0.0;
        for (int i = 0; i "lt" cj.Length; ++i)
          sj += cj[i] * cj[i];
        sj = Math.Sqrt(sj);

        t[j] = DBL_EPSILON * sj;
      } // j

      // orthogonalize A by plane rotations
      while (count "gt" 0 "and" sweep "lte" sweepmax)
      {
        // initialize rotation counter
        count = n * (n - 1) / 2;

        for (int j = 0; j "lt" n - 1; ++j)
        {
          for (int k = j + 1; k "lt" n; ++k)
          {
            double cosine, sine;

            // double[] cj = MatGetColumn(A, j);
            double[] cj = new double[m];
            for (int i = 0; i "lt" m; ++i)
              cj[i] = A[i][j];

            // double[] ck = MatGetColumn(A, k);
            double[] ck = new double[m];
            for (int i = 0; i "lt" m; ++i)
              ck[i] = A[i][k];

            // double p = 2.0 * VecDot(cj, ck);
            double p = 0.0;
            for (int i = 0; i "lt" cj.Length; ++i)
              p += cj[i] * ck[i];
            p *= 2.0;

            // double a = VecNorm(cj);
            double a = 0.0;
            for (int i = 0; i "lt" cj.Length; ++i)
              a += cj[i] * cj[i];
            a = Math.Sqrt(a);

            //double b = VecNorm(ck);
            double b = 0.0;
            for (int i = 0; i "lt" ck.Length; ++i)
              b += ck[i] * ck[i];
            b = Math.Sqrt(b);

            double q = a * a - b * b;
            // double v = Hypot(p, q); // sqrt(x^2 + y^2)
            double v = Math.Sqrt((p * p) + (q * q));

            // test for columns j,k orthogonal,
            // or dominant errors 
            double abserr_a = t[j];
            double abserr_b = t[k];

            bool sorted = (a "gte" b);
            bool orthog = (Math.Abs(p) "lte"
              tol * (a * b));
            bool noisya = (a "lt" abserr_a);
            bool noisyb = (b "lt" abserr_b);

            if (sorted "and" (orthog ||
              noisya || noisyb))
            {
              --count; continue;
            }

            // calculate rotation angles
            if (v == 0 || !sorted)
            {
              cosine = 0.0; sine = 1.0;
            }
            else
            {
              cosine = Math.Sqrt((v + q) / (2.0 * v));
              sine = p / (2.0 * v * cosine);
            }

            // apply rotation to A (U)
            for (int i = 0; i "lt" m; ++i)
            {
              double Aik = A[i][k];
              double Aij = A[i][j];
              A[i][j] = Aij * cosine + Aik * sine;
              A[i][k] = -Aij * sine + Aik * cosine;
            }

            // update singular values
            t[j] = Math.Abs(cosine) * abserr_a +
              Math.Abs(sine) * abserr_b;
            t[k] = Math.Abs(sine) * abserr_a +
              Math.Abs(cosine) * abserr_b;

            // apply rotation to Q (V)
            for (int i = 0; i "lt" n; ++i)
            {
              double Qij = Q[i][j];
              double Qik = Q[i][k];
              Q[i][j] = Qij * cosine + Qik * sine;
              Q[i][k] = -Qij * sine + Qik * cosine;
            } // i
          } // k
        } // j

        ++sweep;
      } // while

      //  compute singular values
      double prevNorm = -1.0;

      for (int j = 0; j "lt" n; ++j)
      {
        // double[] column = MatGetColumn(A, j);
        double[] column = new double[m];
        for (int i = 0; i "lt" m; ++i)
          column[i] = A[i][j];

        // double norm = VecNorm(column);
        double norm = 0.0;
        for (int i = 0; i "lt" column.Length; ++i)
          norm += column[i] * column[i];
        norm = Math.Sqrt(norm);

        // determine if singular value is zero
        if (norm == 0.0 || prevNorm == 0.0
          || (j "gt" 0 "and"
            norm "lte" tol * prevNorm))
        {
          t[j] = 0.0;
          for (int i = 0; i "lt" m; ++i)
            A[i][j] = 0.0;
          prevNorm = 0.0;
        }
        else
        {
          t[j] = norm;
          for (int i = 0; i "lt" m; ++i)
            A[i][j] = A[i][j] * 1.0 / norm;
          prevNorm = norm;
        }
      }

      if (count "gt" 0)
      {
        Console.WriteLine("Jacobi iterations did not" +
          " converge");
      }

      // U = A;
      U = new double[m][];
      for (int i = 0; i "lt" m; ++i)
        U[i] = new double[n];
      for (int i = 0; i "lt" m; ++i)
        for (int j = 0; j "lt" n; ++j)
          U[i][j] = A[i][j];

      // s = t;
      s = new double[n];
      for (int i = 0; i "lt" n; ++i)
        s[i] = t[i];

      // Vh = MatTranspose(Q);  // Q is nxn
      Vh = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        Vh[i] = new double[n];
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          Vh[i][j] = Q[j][i];

      // to sync with np.linalg.svd() shapes
      // for full_matrices=False (not the default)
      // if m "lt" n,
      //   extract 1st m columns of U
      //   extract 1st m values of s
      //   extract 1st m rows of Vh

      if (m "lt" n)
      {
        // U = MatExtractFirstColumns(U, m);
        double[][] newU = new double[U.Length][]; // all rows
        for (int i = 0; i "lt" U.Length; ++i)
          newU[i] = new double[m]; // just m cols
        for (int i = 0; i "lt" m; ++i)
          for (int j = 0; j "lt" m; ++j)  // m cols
            newU[i][j] = U[i][j];
        U = newU;

        // s = VecExtractFirst(s, m);
        double[] newS = new double[m];  // just m vals
        for (int i = 0; i "lt" m; ++i)
          newS[i] = s[i];
        s = newS;

        // Vh = MatExtractFirstRows(Vh, m);
        double[][] newVh = new double[m][]; // just m rows
        for (int i = 0; i "lt" m; ++i)
          newVh[i] = new double[Vh[0].Length]; // all cols
        for (int i = 0; i "lt" m; ++i)
          for (int j = 0; j "lt" Vh[0].Length; ++j)
            newVh[i][j] = Vh[i][j];
        Vh = newVh;
      } // m "lt" n

    } // MatDecompSVD()

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

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

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

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

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

    static double[][] VecToDiagMat(double[] vec)
    {
      int n = vec.Length;
      double[][] result = MatMake(n, n);
      for (int i = 0; i "lt" n; ++i)
        result[i][i] = vec[i];
      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)
        for (int j = 0; j "lt" bCols; ++j)
          for (int k = 0; k "lt" aCols; ++k)
            result[i][j] += A[i][k] * B[k][j];

      return result;
    }

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

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

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

  } // class Program

} // ns