If you have a square n-by-n matrix A, the matrix has a characteristic polynomial with n+1 coefficients (the first coefficient is always 1).

For example, if A =

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

then the coefficients of the characteristic polymonial of A are (1.0, -15.0, -4.0, 231.0, -192.0, -385.0).

The most common technique to compute the characteristic coefficients is to use the Faddeev LeVerrier method. However, the FV method is very unstable in the sense it will blow up for even small matrices that have all positive values because the method computes A*A, A*A*A, A*A*A*A, etc.

An alternative technique to compute the characteristic coefficients is to use polynomial expansion based on the eigenvalues of A. For example, if some matrix A has eigenvalues (2, -3, 4, 5) you would set p(x) = (x-2) * (x+3) * (x-4) * (x-5) and expand using FIFO to get

p(x) = (x-2)*(x+3)*(x-4)*(x-5)
     = x^4 - 8x^3 + 5x^2 + 74x - 120

and this gives you the coefficients = (1, -8, 5, 74, -120. The polynomial exapansion technique has the advantage of being stable, but the significant disadvantage of requiring the eigenvalues.

A few days ago, I implemented a function to compute the coefficients of the characteristic polynomial of a matrix using polynomial expansion, using the Python language. One Sunday morning, I decided to refactor that code to C#.

I quickly ran into a problem I hadn’t considered. The eigenvalues of a matrix can be real or complex (involving the square root of -1). Python handles complex numbers automatically but C# does not. So, to simplify my exploration using C#, I made the assumption that the source matrix is symmetric. The eigenvalues of a symmetric matrix are all real.

With that simplifying assumption, my demo program had no problems. The most difficult part of the demo program is computing the eigenvalues of the source matrix. I used the QR technique, but I could have used the Jacobi technique.

The output of the demo is:

Begin coefficients of characteristic polynomial using
  polynomial expansion.
Square symmetric matrices only.

Source matrix A:
   1.0   2.0   3.0   1.0   5.0
   2.0   5.0   4.0   1.0   4.0
   3.0   4.0   0.0   2.0   2.0
   1.0   1.0   2.0   3.0   1.0
   5.0   4.0   2.0   1.0   6.0

Eigenvalues of A:
  14.2144  -3.7908   3.1914   2.3412  -0.9563

Coefficents of characteristic polynomial:
    1.0  -15.0  -4.0  231.0  -192.0  -385.0

End demo

I validated my demo code using the NumPy library built-in np.poly() function.

It was an interesting and fun way to start a Sunday. Not a particularly useful algorithm, but interesting.

“The Adventures of Rocky and Bullwinkle” (1959-1964) was an animated series. In addition to episodes starring Rocky (the flying squirrel) and Bullwinkle (the moose), there were other segments. All were aimed at both children and adults. Their primary characteristic is that they relied on sophisticated dialog rather than sophisticated art. I enjoyed all the segments of the Rocky and Bullwinkle show as a young man, and I enjoy them even today.

“Fractured Fairy Tales” featured retellings of classic fairy tales with altered storylines, often modernized or humorous, to present a new perspective or message.

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

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

// compute the coefficients of the characteristic
// polynomial of a square matrix using polynomial
// expansion instead of usual Faddeev-LeVerrier algorithm

// for symmetric matrices only - all eigenvalues are 
// real, no complex.

// 1. compute eigenvalues, 2.) compute coefficients

namespace MatCharacteristicPolyExpansion
{
  internal class MatCharacteristicPolyExpansionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin coefficients of " +
        "characteristic polynomial using polynomial " +
        "expansion ");
      Console.WriteLine("Square symmetric matrices only.");

      double[][] A = new double[5][];
      A[0] = new double[] { 1.0, 2.0, 3.0, 1.0, 5.0 };
      A[1] = new double[] { 2.0, 5.0, 4.0, 1.0, 4.0 };
      A[2] = new double[] { 3.0, 4.0, 0.0, 2.0, 2.0 };
      A[3] = new double[] { 1.0, 1.0, 2.0, 3.0, 1.0 };
      A[4] = new double[] { 5.0, 4.0, 2.0, 1.0, 6.0 };

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

      double[] eigenvals = Utils.Eigenvalues(A);
      Console.WriteLine("\nEigenvalues of A: ");
      Utils.VecShow(eigenvals, 4, 9);

      double[] coeffs = Utils.MatPoly(A);
      Console.WriteLine("\nCoefficents of characteristic " +
        "polynomial: ");
      Utils.VecShow(coeffs, 1, 7);

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

  } // class Program

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

  public class Utils
  {
    // ------------------------------------------------------

    public static double[] PolyFromEigenvalues(double[] evals)
    {
      // coefficients of characteristic polynomial
      // from eigenvalues
      // eigenvalues must all be real -- symmetric square src
      int n = evals.Length;
      double[] old = new double[n + 1];
      double[] next = new double[n + 1];

      old[0] = 1.0;
      old[1] = -1 * (evals[0] + evals[1]);
      old[2] = evals[0] * evals[1];

      for (int i = 0; i "lt" n - 2; ++i)
      {
        double currEigen = evals[i + 2];
        next[0] = 1.0;
        for (int j = 1; j "lt" i + 3; ++j)
          next[j] = old[j] + (-1 * old[j - 1] * currEigen);
        next[i + 3] = old[i + 2] * (-1 * currEigen);
        for (int k = 0; k "lt" n + 1; ++k)
          old[k] = next[k];
      }
      return next;
    }

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

    public static double[] MatPoly(double[][] A)
    {
      double[] eigenvals = Eigenvalues(A);
      double[] result = PolyFromEigenvalues(eigenvals);
      return result;
    }

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

    public static double[] Eigenvalues(double[][] A)
    {
      // eigenvalues of square symmetric matix
      // using QR decomposition
      int n = A.Length;
      double[][] X = MatCopyOf(A);  // square, symmetric
      double[][] Q; double[][] R;
      int maxCt = 10000;

      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          if (A[i][j] != A[j][i])
            Console.WriteLine("FAIL. Matrix must" +
              " be square symmetric ");

      int ct = 0;
      while (ct "lt" maxCt)
      {
        MatDecomposeQR(X, out Q, out R, false);
        X = MatMult(R, Q);  // note order
        ++ct;

        if (MatIsUpperTri(X, 1.0e-8) == true)
          break;
      }

      if (ct "gte" maxCt)
        Console.WriteLine("Warn: Eigenvalues convergence ");

      double[] result = new double[n];
      for (int j = 0; j "lt" n; ++j)
        result[j] = X[j][j];
      return result;
    }

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

    public static void MatDecomposeQR(double[][] mat,
      out double[][] q, out double[][] r,
      bool standardize)
    {
      // QR decomposition, Householder algorithm.
      // assumes square matrix

      int n = mat.Length;  // assumes mat is nxn
      int nCols = mat[0].Length;
      if (n != nCols) Console.WriteLine("M not square ");

      double[][] Q = MatIdentity(n);
      double[][] R = MatCopyOf(mat);
      for (int i = 0; i "lt" n - 1; ++i)
      {
        double[][] H = MatIdentity(n);
        double[] a = new double[n - i];
        int k = 0;
        for (int ii = i; ii "lt" n; ++ii)  // last part col [i]
          a[k++] = R[ii][i];

        double normA = VecNorm(a);
        if (a[0] "lt" 0.0) { normA = -normA; }
        double[] v = new double[a.Length];
        for (int j = 0; j "lt" v.Length; ++j)
          v[j] = a[j] / (a[0] + normA);
        v[0] = 1.0;

        double[][] h = MatIdentity(a.Length);
        double vvDot = VecDot(v, v);
        double[][] alpha = VecToMat(v, v.Length, 1);
        double[][] beta = VecToMat(v, 1, v.Length);
        double[][] aMultB = MatMult(alpha, beta);

        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            h[ii][jj] -= (2.0 / vvDot) * aMultB[ii][jj];

        // copy h into lower right of H
        int d = n - h.Length;
        for (int ii = 0; ii "lt" h.Length; ++ii)
          for (int jj = 0; jj "lt" h[0].Length; ++jj)
            H[ii + d][jj + d] = h[ii][jj];

        Q = MatMult(Q, H);
        R = MatMult(H, R);
      } // i

      if (standardize == true)
      {
        // standardize so R diagonal is all positive
        double[][] D = MatMake(n, n);
        for (int i = 0; i "lt" n; ++i)
        {
          if (R[i][i] "lt" 0.0) D[i][i] = -1.0;
          else D[i][i] = 1.0;
        }
        Q = MatMult(Q, D);
        R = MatMult(D, R);
      }

      q = Q;
      r = R;

    } // QR decomposition

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

    public static bool MatIsUpperTri(double[][] mat,
      double tol)
    {
      int n = mat.Length;
      for (int i = 0; i "lt" n; ++i)
      {
        for (int j = 0; j "lt" i; ++j)
        {  // check lower vals
          if (Math.Abs(mat[i][j]) "gt" tol)
          {
            return false;
          }
        }
      }
      return true;
    }

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

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

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

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

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

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

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

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

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

    public static double[][] MatAdd(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;
    }

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

    public static double MatTrace(double[][] A)
    {
      double result = 0.0;
      for (int i = 0; i "lt" A.Length; ++i)
        result += A[i][i];
      return result;
    }

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

    public static double[][] MatCopyOf(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;
    }

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

    public static double[][] MatScalarMult(double u,
      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] = u * A[i][j];
      return result;
    }

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

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

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

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

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

    public static double[][] MatRandom(int dim, Random rnd)
    {
      // dim-by-dim matrix with random values in [-1.0, 1.0)
      double[][] result = MatMake(dim, dim);
      double lo = -1.0; double hi = 1.0;
      for (int i = 0; i "lt" dim; ++i)
        for (int j = 0; j "lt" dim; ++j)
          result[i][j] = (hi - lo) * rnd.NextDouble() + lo;
      return result;
    }

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

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

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

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

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

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

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

    public 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-8) x = 0.0;
        Console.Write(x.ToString("F" + dec).
          PadLeft(wid));
      }
      Console.WriteLine("");
    }

  } // class Utils

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

} // ns