In mathematics, a condition number is a measure of how sensitive an object (vector, array, etc.) is to some operation (norm, inverse, etc.) There are many kinds of condition numbers. The condition number I come across most often is the condition number of a matrix for the inverse operation.

If a matrix has a small condition number, it’s not sensitive, and most inverse algorithms will succeed. If a condition number if large, it is sensitive and fragile algorithms may fail.

The condition number for matrix inverse is usually (but not always) defined as the ratio of the largest singular value to the smallest singular value.

So, the straightforward approach to compute a matrix condition number for matrix inverse is to compute the singular values using one of several singular value decomposition algorithms (the Jacobi technique and the QR technique are common), then find the largest and smallest singular values and divide.

A few days ago, I did an experiment to compute the matrix inverse condition number in a very convoluted way. I used Power Iteration and Rayleigh Quotient to compute the largest eigenvalue and then the largest singular value (square root of the eigenvalue) of the source matrix. Then I used Newton Iteration inverse to find the inverse of the source matrix, and then Power Iteration and Rayleigh Quotient on the inverse to find the smallest singular value of the source matrix. It was an interesting exploration, but not very practical for several reasons. For example, eigenvalues can be complex numbers (involving the square root of -1) in which case the whole technique becomes a mess.

I was sitting in a hotel lobby in Boston, waiting for my room to get ready, and figured I’d write some code to compute the condition number for matrix inverse, using SVD (Jacobi version) directly.

A couple of notes before I forget. Because I’m computing a condition number for matrix inverse, the source matrix must be square. The source matrix does not have to be symmetric — there are specialized algorithms to compute the singular value decomposition of symmetric matrices. Singular values of a square matrix are always non-negative. The smallest singular value of a square matrix could be zero but only when the determinant of the source matrix is zero, which means the source matrix has no inverse.

Anyway, I had a demo up and running pretty quickly. The key calling statements are:

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

double[][] U;
double[][] Vh;
double[] s;
Utils.SVD_Jacobi(A, out U, out Vh, out s);

double largeSingular = Utils.VecMax(s);
double smallSingular = Utils.VecMin(s);
double cn = largeSingular / smallSingular; // no 0.0

Console.WriteLine("\nCondition number of A = " +
  cn.ToString("F4"));

The output of the demo is:

Begin matrix condition number using SVD-Jacobi

Source matrix A:
   3.0000   0.0000  -2.0000   5.0000
  -1.0000   4.0000   6.0000   3.0000
   4.0000   1.0000   0.0000   3.0000
  -3.0000   2.0000   4.0000   5.0000

Calling SVD()

Singular values:
  10.5124   7.5957   3.4206   0.3075

Large singular = 10.5124
Small singular = 0.3075

Condition number of A = 34.1821

End demo

The condition number, 34.1821, is relatively small, so the source matrix isn’t very sensitive to matrix inverse, and most inverse algorithms should work.

In my opinion, condition numbers are more of an interesting theory idea than they are practical techniques. But I live in the world of machine learning, and condition numbers might be important and practical in fields I’m not familiar with.

Much fun, and as I’m writing this sentence, my hotel room at the Copley Square Hotel just got ready. Win-win.

I love baseball. I have fond memories of my father taking my brother and me to see the California/Anaheim Angels just a few miles away from our house in neighboring Fullerton.

Right: I’ve never been to Boston’s Fenway Park. But it might be worth it just to see the infamous worst seat in Major League Baseball. I would give this seat a very poor condition number.

Demo code. 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;

// matrix condition number using SVD - Jacobi version
// find singular values, ratio of largest to smallest

namespace MatrixConditionNumberSVDJacobi
{
  internal class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin matrix condition number" +
        " using SVD-Jacobi ");

      // 1. set up source matrix, must be square
      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 };

      Console.WriteLine("\nSource matrix A: ");
      Utils.MatShow(A, 4, 9);

      double[][] U;
      double[][] Vh;
      double[] s;
      Utils.SVD_Jacobi(A, out U, out Vh, out s);
      Console.WriteLine("\nCalling SVD() ");
      Console.WriteLine("\nSingular values: ");
      Utils.VecShow(s, 4, 9);

      double largeSingular = Utils.VecMax(s);
      double smallSingular = Utils.VecMin(s);
      Console.WriteLine("\nLarge singular = " +
        largeSingular.ToString("F4"));
      Console.WriteLine("Small singular = " +
        smallSingular.ToString("F4"));
      double cn = largeSingular / smallSingular; // no 0.0

      Console.WriteLine("\nCondition number of A = " +
        cn.ToString("F4"));

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

  } // class Program

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

    public static double MatConditionNumber(double[][] M)
    {
      // assumes M is square
      double[][] U; double[][] Vh; double[] s;
      SVD_Jacobi(M, out U, out Vh, out s);
      double large = VecMax(s);
      double small = VecMin(s);  // assume not 0
      return large / small;
    }

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

    public static void SVD_Jacobi(double[][] M,
      out double[][] U, out double[][] Vh, out double[] s)
    {
      // github.com/ampl/gsl/blob/master/linalg/svd.c
      double DBL_EPSILON = 1.0e-15;

      double[][] A = MatCopyOf(M); // working U
      int m = A.Length;
      int n = A[0].Length;
      double[][] Q = MatIdentity(n); // working V
      double[] t = new double[n];  // working s

      // initialize counters
      int count = 1;
      int sweep = 0;
      //int sweepmax = 5 * n;

      double tolerance = 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 = Utils.MatGetColumn(A, j);
        double sj = Utils.VecNorm(cj);
        t[j] = DBL_EPSILON * sj;
      }

      // orthogonalize A by plane rotations
      while (count "gt" 0 && 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[] ck = MatGetColumn(A, k);

            double p = 2.0 * VecDot(cj, ck);
            double a = VecNorm(cj);
            double b = VecNorm(ck);

            double q = a * a - b * b;
            double v = Hypot(p, 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"
              tolerance * (a * b));
            bool noisya = (a "lt" abserr_a);
            bool noisyb = (b "lt" abserr_b);

            if (sorted && (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 norm = VecNorm(column);

        // determine if singular value is zero
        if (norm == 0.0 || prevNorm == 0.0
          || (j "gt" 0 &&
            norm "lte" tolerance * 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;
      Vh = MatTranspose(Q);
      s = t;

      // to sync with default np.linalg.svd() shapes:
      // 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);
        s = VecExtractFirst(s, m);
        Vh = MatExtractFirstRows(Vh, m);
      }

    } // SVD_Jacobi()

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

    public static double[] MatGetColumn(double[][] M, int j)
    {
      int nRows = M.Length;
      double[] result = new double[nRows];
      for (int i = 0; i "lt" nRows; ++i)
        result[i] = M[i][j];
      return result;
    }

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

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

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

    public static double Hypot(double x, double y)
    {
      // fancy sqrt(x^2 + y^2)
      double xabs = Math.Abs(x);
      double yabs = Math.Abs(y);
      double min, max;

      if (xabs "lt" yabs)
      {
        min = xabs; max = yabs;
      }
      else
      {
        min = yabs; max = xabs;
      }

      if (min == 0)
        return max;

      double u = min / max;
      return max * Math.Sqrt(1 + u * u);
    }

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

    public static double[][] MatExtractFirstColumns(double[][] M,
      int n)
    {
      int nRows = M.Length;
      // int nCols = src[0].Length;
      double[][] result = MatMake(nRows, n);
      for (int i = 0; i "lt" nRows; ++i)
        for (int j = 0; j "lt" n; ++j)
          result[i][j] = M[i][j];
      return result;
    }

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

    public static double[][] MatExtractFirstRows(double[][] M,
      int n)
    {
      // int nRows = src.Length;
      int nCols = M[0].Length;
      double[][] result = MatMake(n, nCols);
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" nCols; ++j)
          result[i][j] = M[i][j];
      return result;
    }

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

    public static double[] VecExtractFirst(double[] v, int n)
    {
      double[] result = new double[n];
      for (int i = 0; i "lt" n; ++i)
        result[i] = v[i];
      return result;
    }

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

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

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

    public static double VecMaxMagnitude(double[] v)
    {
      double mx = Math.Abs(v[0]);
      for (int i = 0; i "lt" v.Length; ++i)
        if (Math.Abs(v[i]) "gt" mx)
          mx = v[i];
      return mx;
    }

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

    public static double VecMax(double[] v)
    {
      double mx = v[0];
      for (int i = 0; i "lt" v.Length; ++i)
        if (v[i] "gt" mx)
          mx = v[i];
      return mx;
    }

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

    public static double VecMin(double[] v)
    {
      double mn = v[0];
      for (int i = 0; i "lt" v.Length; ++i)
        if (v[i] "lt" mn)
          mn = v[i];
      return mn;
    }

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

    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 void VecShow(double[] v,
      int dec, int wid)
    {
      for (int i = 0; i "lt" v.Length; ++i)
      {
        double x = v[i];
        if (Math.Abs(x) "lt" 1.0e-8) x = 0.0;
        Console.Write(x.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[][] 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[][] MatTranspose(double[][] A)
    {
      int nr = A.Length; int nc = A[0].Length;
      double[][] result = MatMake(nc, nr);  // note reverse
      for (int i = 0; i "lt" nr; ++i)
        for (int j = 0; j "lt" nc; ++j)
          result[j][i] = 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;
    }

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

  } // class Utils
} // ns