I’ve been working with machine learning using the C# language recently. A fairly common operation is multiplying a matrix and a vector together. The topic is surprisingly tricky.

First, I created and ran a demo using Python/NumPy. If A is a matrix and v is a vector, matmul(A, v) and matmul(v, A) both return a vector even though the return could be a single-column matrix or a single-row matrix.

Next, I implemented MatVecProd(A, v) and VecMatProd(v, A) functions using C#. I’m a good coder, and I have a lot of experience working with C# vectors and matrices, but even so, I made a lot of mistakes. The effort I thought would take about 20 minutes ended up taking a couple of hours.

An alternative strategy for multiplying a matrix and a vector together is to convert the vector to a matrix (with either one row or one column, depending on the scenario), multiply, and then either use the resulting matrix (with one row or one column) or converting the result matrix back to a vector. I often use this technique because it’s more clear and less error-prone, but the technique is much less efficient than a direct approach.

There’s no major moral to the story here, but I guess the minor moral is that sometimes implementing an algorithm is more difficult than expected.

When multiplying a matrix and a vector, I was a bit surprised at how difficult it was to combine two different but similar data structures. In science fiction movies, there are quite a few examples that combine a human and a robot.

Right: In “Virus” (1999), the crew of a salvage ship are hit by a storm and seek refuge on a Soviet research ship. Unfortunately, the research ship has been taken over by an alien entity of some sort and it sets out to capture and robotize the crew. This movie was a huge flop but I inexplicably like it quite a bit. I give it a B+ grade but most critics gave it a D grade.

Demo program. Replace “lt” with Boolean less-than operator symbol.

using System;
using System.IO;

namespace MatrixLibrary
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin matrix * vector demo ");

      double[][] A = MatCreate(2, 3);
      A[0] = new double[] { 1, 2, 3 };
      A[1] = new double[] { 4, 5, 6 };
      Console.WriteLine("\nmatrix A: ");
      MatShow(A, 0, 5);

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

      double[] v = new double[] { 7, 8, 9 };
      Console.WriteLine("\nvector v: ");
      VecShow(v, 0, 5);

      double[] Av = MatVecProd(A, v);
      Console.WriteLine("\nA * v as vector = ");
      VecShow(Av, 0, 5);

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

      double[] w = new double[] { 7, 8 };
      Console.WriteLine("\nvector w: ");
      VecShow(w, 0, 5);

      double[] wA = VecMatProd(w, A);
      Console.WriteLine("\nw * A as vector = ");
      VecShow(wA, 0, 5);

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

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

    public static double[] MatVecProd(double[][] m,
      double[] v)
    {
      // return a regular vector
      int nRows = m.Length; int nCols = m[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] += m[i][k] * v[k];

      return result;
    }

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

    public static double[] VecMatProd(double[] v,
      double[][] m)
    {
      // return a regular vector
      int nRows = m.Length; int nCols = m[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] * m[k][j];

      return result;
    }

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

    static double[][] MatCreate(int rows, int cols)
    {
      double[][] result = new double[rows][];
      for (int i = 0; i "lt" rows; ++i)
        result[i] = new double[cols];
      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];
          if (Math.Abs(v) "lt" 1.0e-5) v = 0.0;  // hack
          Console.Write(v.ToString("F" + dec).PadLeft(wid));
        }
        Console.WriteLine("");
      }
    }

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

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

  } // Program
} // ns