In mathematics, a condition number is a measure of how sensitive an object (array, matrix, etc.) is to some operation (dot product, matrix inverse, etc.) There are many different kinds of condition numbers. The condition number I come across most often is one that computes sensitivity of matrix inverse.
The condition number version that I use is simple: the ratio of the largest singular value of the matrix to the smallest. For example, if 4-by-4 matrix A is:
[[ 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]]
The four singular values are (10.5124, 7.5957, 3.4206, 0.3075). The largest singular value is 10.5124 and the smallest is 0.3075 and so the condition number is 10.5124 / 0.3075 = 34.1821. This is a relatively small value so the matrix is not very sensitive to the inverse operation and most inverse algorithms should work.
So, to compute a condition number for matrix inverse sensitivity, all you need to do is compute the singular values of the source matrix. There are several different algorithms, but for some subjective reason, on this morning, I just didn’t want to use an SVD algorithm.
So, one morning before work, I made an inexplicable decision to compute the condition number for matrix inverse using eigenvalues. Well, computing eigenvalues is even more difficult and ugly than computing singular values. However, there is a shortcut algorithm to compute just the largest and smallest eigenvalues. If you have a source matrix A, you compute At * A where At is the transpose of A (easy). Then you use Power Iteration and Rayleigh quotient to compute the largest eigenvalue of At * A (relatively easy). Then you use Newton iteration to find the inverse of At * A (relatively easy), followed by Power Iteration and Rayleigh quotient to find the smallest eigenvalue of At * A (relatively easy). Now if you take the square roots of the largest and smallest eigenvalue of At * A (super easy), you get the largest and smallest singular values of the source matrix A. Then you divide the larger singular value by the smaller singular value to get the condition number for matrix inverse of A.
Put another way, the steps above compute the largest and smallest singular values of a source matrix in a series of relatively easy steps rather than computing the singular values directly using SVD which is difficult.
I implemented a demo in C#. The output is:
Begin matrix condition number 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 At * A: 35.0000 -6.0000 -24.0000 9.0000 -6.0000 21.0000 32.0000 25.0000 -24.0000 32.0000 56.0000 28.0000 9.0000 25.0000 28.0000 68.0000 Largest eigenvalue of At * A = 110.5096 Smallest eigenvalue of At * A = 0.0946 Singular values of A = 10.5124, 0.3075 Condition number of A for matrix inverse = 34.1821 End demo
I verified the results using a Python language program with the scipy.linalg.svdvals() function that computes singular values. Let me emphasize this was a thought experiment, not a practical technique.
A very interesting and challenging exploration to start the day.
On the night before I wrote this blog post, I was at a friend’s house with a handful of people. We played this horse race dice game. It was a lot of fun! Notice that Horse #2 only needs a few rolls of the dice to win but Horse #7 needs many rolls because 7 is so common.
Rules:
From the deck of cards, discard the aces, kings, and jokers, leaving you with 44 cards. Shuffle and deal these cards. Each player will not necessarily receive the same amount of cards.
Place the horses at the starting gate. Four horses must be eliminated, or Scratched before the race starts. To do this, the player to the left of the dealer rolls the dice, and the total on the dice is the first horse to be scratched. That horse is moved back to the first line. For example, if a three and a five were thrown with the dice, a total of eight, then horse number eight would be moved back to the first line. All players would then check their cards and pay five cents for each number eight card in their hand. The eights are then discarded.
The next player then rolls the dice for the second scratch. That horse moves back to the second line and costs ten cents for each of those that you have in your hand. Continue with the third and fourth scratch, moving to the fifteen and twenty cent lines, respectively. Paying the designated amounts, and discarding the scratched cards accordingly. If you shake the number of a horse that is already scratched, you pay the amount on that horse’s line, and pass the dice to the next player.
The horses now remaining at the gate are ready for the race to start. The next person to roll the dice totals them, and moves that horse forward one space. However, if the total on the dice is that of a scratched horse, then that player must pay the amount which that horse is standing on. Play continues around the table, until one horse reaches the finish line!
Your cards are your shares of each horse in the race. If horse number nine should win the race, and you have no nines in your hand, then you’re just plain out of luck. On the other hand, if you were to have two nines in your hand, you would win half the pot!
Each card with the winning horse’s number wins 25% of the pot. If you had all four of the one number, and that horse was to win then you would take the entire pot. After the pot has been divided, the cards are passed to the next player, the horses are all moved to the starting gate, and you’re ready to begin the next race!
Poker chips can be substituted for money. A pair of ones on the dice would be a #2 horse. A pair of sixes would be a Queen. A five and a six would be a Jack.
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;
// matrix condition number using power iteration, Rayleigh
// quotient, and Newton iteration inverse
namespace MatrixConditionNumber
{
internal class MatrixConditionNumberProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nBegin matrix condition number ");
// 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 };
Console.WriteLine("\nSource matrix A: ");
Utils.MatShow(A, 4, 9);
// 2. compute At * A
double[][] At = Utils.MatTranspose(A);
double[][] AtA = Utils.MatMult(At, A);
Console.WriteLine("\nAt * A: ");
Utils.MatShow(AtA, 4, 9);
// 3. compute largest eigenvalue of AtA using PI and RQ
double domEigenvalue =
Utils.MatDomEigenvalue(AtA);
Console.WriteLine("\nLargest eigenvalue of At * A = " +
domEigenvalue.ToString("F4"));
// 4. compute smallest eigenvalue using NI, PI, RQ
double leastDomEigenvalue =
Utils.MatLeastDomEigenvalue(AtA);
Console.WriteLine("Smallest eigenvalue of At * A = " +
leastDomEigenvalue.ToString("F4"));
// 5. compute condition number
double largeSV = Math.Sqrt(domEigenvalue);
double smallSV = Math.Sqrt(leastDomEigenvalue);
Console.WriteLine("\nSingular values of A = " +
largeSV.ToString("F4") + ", " +
smallSV.ToString("F4"));
double conditionNumber = largeSV / smallSV;
Console.WriteLine("\nCondition number of A for" +
" matrix inverse = " +
conditionNumber.ToString("F4"));
Console.WriteLine("\nEnd demo ");
Console.ReadLine();
} // Main
} // class Program
// ========================================================
public class Utils
{
// ------------------------------------------------------
public static double MatDomEigenvalue(double[][] A,
int maxIter = 10000, double tol = 1.0e-8)
{
double[] domEigenvector =
MatDomEigenvector(A, maxIter, tol);
double domEigenvalue =
MatEigenvalueFromEigenvector(A, domEigenvector);
return domEigenvalue;
}
// ------------------------------------------------------
public static double MatLeastDomEigenvalue(double[][] A,
int maxIter = 10000, double tol = 1.0e-8)
{
double[][] Ainv =
MatInverseNewton(A, maxIter, tol);
// check inverse
double[][] C = MatMult(A, Ainv);
double[][] I = MatIdentity(A.Length);
//if (MatAreEqual(C, I, 1.0e-8) == true)
// Console.WriteLine("Newton inverse OK ");
//else
// Console.WriteLine("Newton inverse FAIL ");
double domEigenvalue =
MatDomEigenvalue(Ainv, maxIter, tol);
double leastDomEigenvalue = 1.0 / domEigenvalue;
return leastDomEigenvalue;
}
// ------------------------------------------------------
public static double[] MatDomEigenvector(double[][] A,
int maxIter, double tol)
{
// power iteration technique
int n = A.Length; ;
double[] xCurr = new double[n];
for (int i = 0; i "lt" n; ++i)
xCurr[i] = 1.0;
double[] xNew = VecScalarProduct(2.0, xCurr);
int iter = 0;
while (iter "lt" maxIter)
{
xNew = Utils.MatVecProd(A, xCurr);
double scale = Utils.VecMaxMagnitude(xNew);
xNew = Utils.VecScalarProduct(1.0 / scale, xNew);
if (Utils.VecAreEqual(xCurr, xNew, tol) == true)
break;
xCurr = xNew;
++iter;
}
// Console.WriteLine("MatDominantEigenvector needed " +
// iter + " iterations ");
if (iter "gte" maxIter)
Console.WriteLine("WARN: no converge ");
return xCurr;
}
// ------------------------------------------------------
public static double
MatEigenvalueFromEigenvector(double[][] A,
double[] eigVec)
{
// Rayleigh quotient technique
// (Ax * x) / (x * x)
double[] x = eigVec; // to match documentation
double[] Ax = MatVecProd(A, x);
double Axx = VecDot(Ax, x);
double xx = VecDot(x, x);
return Axx / xx; // asumews xx not 0
}
// ------------------------------------------------------
public static double[] VecScalarProduct(double u,
double[] v)
{
int n = v.Length;
double[] result = new double[n];
for (int i = 0; i "lt" n; ++i)
result[i] = v[i] * u;
return result;
}
// ------------------------------------------------------
public static double[] VecMatProd(double[] v,
double[][] A)
{
// one-dim vec * two-dim mat
int nRows = A.Length;
int nCols = A[0].Length;
int n = v.Length;
if (n != nCols)
throw new Exception("non-comform in VecMatProd");
double[] result = new double[n];
for (int i = 0; i "lt" n; ++i)
{
for (int j = 0; j "lt" nCols; ++j)
{
result[i] += v[j] * A[i][j];
}
}
return result;
}
// ------------------------------------------------------
public static bool VecAreEqual(double[] v1, double[] v2,
double tol)
{
int n = v1.Length;
for (int i = 0; i "lt" n; ++i)
{
if (Math.Abs(v1[i] - v2[i]) "gt" tol)
return false;
}
return true;
}
// ------------------------------------------------------
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[] MatVecProd(double[][] A,
double[] v)
{
return VecMatProd(v, A);
}
// ------------------------------------------------------
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[][] 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[][] MatInverseNewton(double[][] A,
int maxIter = 1000, double eps = 1.0e-8)
{
// Newton's method
// X_k+1 = X_k * (2I - A*X_k)
int n = A.Length; // must be square martix
double[][] Xcurr = MatStart(A); // Pan-Reif algorithm
double[][] Xnew = MatMake(n, n);
double[][] I = MatIdentity(n);
double[][] I2 = MatScalarMult(I, 2.0);
int iter = 0;
while (iter "lt" maxIter)
{
Xnew =
MatMult(Xcurr, MatSubtract(I2,
MatMult(A, Xcurr)));
// Xcurr = MatCopyOf(Xnew); // by value
Xcurr = Xnew; // by ref
if (iter % 20 == 0)
{
double[][] check = MatMult(A, Xcurr); //
if (MatAreEqual(check, I, eps) == true)
return Xcurr;
}
++iter;
} // while
Console.WriteLine("WARN: no converge ");
return Xnew;
}
// ------------------------------------------------------
public static double[][] MatStart(double[][] A)
{
// intial matrix for Newton iteration inverse
int n = A.Length;
double maxRowSum = 0.0;
double maxColSum = 0.0;
for (int i = 0; i "lt" n; ++i)
{
double rowSum = 0.0;
for (int j = 0; j "lt" n; ++j)
rowSum += Math.Abs(A[i][j]);
if (rowSum "gt" maxRowSum)
maxRowSum = rowSum;
}
for (int j = 0; j "lt" n; ++j)
{
double colSum = 0.0;
for (int i = 0; i "lt" n; ++i)
colSum += Math.Abs(A[i][j]);
if (colSum "gt" maxColSum)
maxColSum = colSum;
}
double[][] result = MatTranspose(A);
double t = 1.0 / (maxRowSum * maxColSum);
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] *= t;
return result;
}
// ------------------------------------------------------
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[][] MatSubtract(double[][] A,
double[][] B)
{
// A - 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[][] 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[][] A,
double u)
{
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[][] 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;
}
// ------------------------------------------------------
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;
}
} // class Utils
// ========================================================
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2026 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2026 G2E Conference
2026 iSC West Conference
You must be logged in to post a comment.