Ill-conditioned Bases and Numerical Instability

I spent a chunk of this afternoon playing with a MIP model sent to me by someone who could not figure out why a solution CPLEX certified as "optimal" had nontrivial constraint violations.  That motivated today's blog entry, but first a reminiscence (it's my blog, and I grant myself permission).

Back when I started grad school, men were men and computers were mainframes (and had more or less completed the transition from vacuum tubes to transistors; integrated circuits were still pretty experimental).  The first course in my doctoral program was numerical analysis (highly recommended for future OR practitioners), taught by a visiting professor who had not gotten the memo about our "easy grades for doctoral students" policy.  (On the bright side, I spent the remaining 27 quarters of my doctoral program free from worry about maintaining a 4.0 GPA.)  He told an anecdote about a nameless professor who, in the course of some research, inverted a matrix using some homebrewed FORTRAN code. In the inevitable journal article, he printed the coefficients of the inverse matrix (it apparently had a fairly low dimension) to some ridiculous number of decimal places ... all of which were wrong ... as were the first four digits to the left of the decimal point.  The moral of the story was that "stiff" (ill-conditioned) matrices need to be handled carefully.

A bit later in my studies, I took a course in convex analysis that included some linear optimization.  I'd also had a course in linear programming as an undergraduate.  Neither the undergrad course (in which pivoting was done by hand) nor the graduate course (which emphasized theory) mentioned numerical precision at all.  So I'm not surprised that a number of LP/MIP users today are either unaware of or have at best a vague understanding of scaling, conditioning and the perils of stiff matrices.

I probably should offer up some sort of definition of "ill-conditioned" at this point. The Wikipedia has a useful definition of the condition number of a (square, nonsingular) matrix.  Essentially, it's the product of the matrix norm of the matrix with the norm of the inverse matrix. Large condition numbers invite bad rounding errors when dealing with the inverse, and of course the simplex method explicitly or implicitly inverts the basis matrix at each iteration. What constitutes a "large" condition number is in the eye of the beholder, but I think it's safe to say that if your matrix has a condition number bigger than the US federal deficit (in dollars), you should start paying attention to it.

The related concept "numerical instability" generally refers to a situation where small, seemingly benign rounding errors can lead to strange results. One source of small errors is when you export a model to a text file and then import it again. Integers are retrieved at full precision, but decimals are rounded both during export and during import, so the imported model is not identical to the exported model. The differences are small, and if the problem is numerical stable the effect is probably negligible; but if the problem is unstable, those small differences can have significant consequences.

A word about scaling before I move on to stiffness.  If your constraint matrix mixes big honking coefficients and itty-bitty coefficients (I refer here to magnitude, irrespective of sign), there's a pretty good chance that somewhere along the line you'll get some stiff basis matrices.  Even if you don't, there's a good chance rounding error will cause problems somewhere.  To a mathematician, a zillion minus epsilon is less than a zillion (epsilon being small and positive).  To a computer, maybe it's less, maybe it's equal.  So it is advisable to avoid scaling issues by converting the units of variables, where necessary, to bring coefficients into a tighter range of magnitudes.

There is, however, a misconception in some quarters that poor scaling is the only source of ill-conditioned basis matrices. Sadly, that's not true; they can crop up for other, harder to predict, reasons. If you have software that computes condition numbers, pick a nice small value for $\epsilon\gt 0$ and compute the condition number of
\[ \left[\begin{array}{cc} 1 & 1\\ 1-\epsilon & 1 \end{array}\right]. \] Kind of depressing, isn't it?

So what happens if you encounter ill-conditioned basis matrices? It depends on a number of factors, not least of which is the care with which the authors of the software you are using dealt with the possibility. Some algorithms for inverting, factoring or otherwise screwing with matrices are more robust than others, but oftentimes the more robust algorithms are also slower. So it is quite possible (as is the case with CPLEX) that your LP/IP software is capable of compensating for ill-conditioned matrices, at least to some extent, but may not do so to its full capabilities by default because using the full capabilities would slow the software down, including on the majority of problems where that caution was unnecessary.

Optimization software typically will provide, on request, some measures of solution quality. At minimum, you would hope for the maximum violation of any integrality restriction (in a MIP), the maximum violation of any variable bounds, and the maximum violation of any functional constraint (which may be expressed as the maximum bounds violation of any slack variable, which is equivalent). You can, of course, compute these for yourself by substituting the solution values into the constraints. Rounding error being omnipresent, do not expect zero violations across the board; but if you have violations that exceed what you reasonably might expect from rounding, take a look at the condition numbers (sometimes called "kappa values") of the basis matrices. CPLEX will let you see the condition number of the final basis matrix when you solve an LP via one of the simplex variants. A "final" kappa is not available when you solve a MIP, but recent versions of CPLEX allow you, through a parameter setting, to tell it to record the distribution of basis kappa values for either a sample of the node LPs it solves or for all of them (the latter option obviously being slower than the former). It will then tell you what percentage of basis matrices struck it as "stable", "suspect", "unstable" or "ill-posed".

If you are getting suspicious results, there are a number of things you can try, starting with checking the scaling of your coefficients and proceeding to a variety of software settings. I'm familiar with the choices in CPLEX, but my guess is that other optimization software offers a similar palette. First, you can adjust parameter settings that declare the maximum allowable violation of any constraint (EpRHS in CPLEX) or integrality restriction (EpInt). Cranking those tolerances down may help, but it is not guaranteed. (Users sometimes set the tolerances at zero, which can have the unfortunate effect of making a feasible problem look infeasible to CPLEX.) Tighter tolerances may well slow the solution process. Another possibility is to find a setting that tells your solver to try to be more careful about rounding errors. CPLEX has a "Markowitz tolerance" (EpMrk) setting which, if increased, sometimes helps with ill-conditioned bases. The current version (CPLEX 12.2) also has a "numerical precision emphasis" setting (NumericalEmphasis) that tells CPLEX you are concerned about stability issues and want it to be extra careful.