Assume I have a real valued function $f(x_1,\ldots ,x_N)$ of some variables $x_i$ which I want to evaluate numerically. In general the formula for $f$ can contain products, rationals, trancendental functions etc. and will be to long to investigate its numerical stability analytically. Or it will at least be to time consuming to do it in practice. Assume I do not have a shorter equivalent with guaruanteed stability. Is there a methodical procedure to analyse the numerical stability of $f$. I think of comparing it to arbitrary precission results obtained using a computer algebra system. Say the function will be implemented in C using stdlib functions and single or double precision. Which quantities should I compare to quantify the quality of the approximation at finite precission? How do I determine critical values of the variables? How can I choose the compiler and the compiler optimizations so other people can easily reproduce the results? ... I know that the problem setting is probably to generic to give good answers. But I still think that this is a common problem in computer science and wonder if there are references which propose standards to perform such analysis.
3 Answers
What you are looking for is what is called "Automatic error analysis" and is the subject of Chapter 26 of Higham's book "Accuracy and Stability of Numerical Algorithms", 2nd ed., SIAM Publishers.
One technique he describes is using direct search optimization: try to formulate your problem as an optimisation problem and use the optimization algorithm to find coefficients or parameter values that maximize or minimize a quantity related to the accuracy of your algorithm/formula. He uses the example of the growth factor in Gaussian Elimination (what matrix maximizes this growth factor) or the roots of a cubic (as I answered in one of your previous questions).
I would suggest that you obtain a copy of this book, read the introductory chapters and this chapter 26 and the references therein.
Evaluate your function a few times (3 is typically enough) with all inputs slightly randomly perturbed by $\pm 1$ulp. The standard deviation of the three results will give you a crude (but usually sufficient) measure of numerical sensitivity. You can compare this to the expected sensitivity from linearization, and form the quotient to get a stability estimate.
Note that numerical stability asks for how much worse the actual error at evaluation of a particular $x$ is compared to the error expected from sensitivity analysis when changing the inputs by $\pm1$ ulp; the latter error expressed in ulps defines the problem condition. Condition can be very poor for a stable algorithm (example: $1/x$ near $x=0$) and stability can be poor for a very well-conditioned function (example: $1/(1-x)-1/(1+x)$ near $x=0$).
What GertVdE describes is the numerical error. This may be what you are looking for, but it's not the same as numerical stability as indicated in the title of the question. Numerical stability essentially asks whether nearby values of your input variables yield nearby values of the output. In other words, whether a formula like $$|f(x+\varepsilon)-f(x)| \le C |\varepsilon|$$ holds for some moderate constant $C$. For this, you can analyze the derivative of $f$ or, if your function does not vary in its properties wildly over its domain, simply try this for a bunch of $x,\epsilon$ pairs.
-
$\begingroup$ And what can be done if the functions does vary wildly over its domain or if there is no feasible derivative available? Are there other techniques or would we end up with a Monte Carlo approach? $\endgroup$– AndréCommented Sep 7, 2012 at 5:44
-
1$\begingroup$ -1: You explain the notion of condition, not of numerical stability. $\endgroup$ Commented Sep 7, 2012 at 7:17