PARA'04 State-of-the-Art
in Scientific Computing
June 20-23, 2004 (Home page)
Updated: 18 February 2004

Towards joint use of probabilities and intervals in scientific computing: what is the best transition from linear to quadratic approximation?

Martine Ceberio$^1$, Vladik Kreinovich$^1$, and Lev Ginzburg${2,3}$
$^1$ Comp. Sci., U. Texas, El Paso, TX 79968, USA
email: {mceberio,vladik}@cs.utep.edu
$^2$ Dept. of Ecology and Evolution
State University of New York
Stony Brook, NY 11794, USA
email: lev@ramas.com
$^3$ Applied Biomathematics
100 N. Country Road
Setauket, NY 11733, USA

In many areas of science and engineering, we are interested in the value of a physical quantity $y$ that is difficult (or even impossible) to measure directly. Examples may include the amount of a pollutant in a given lake, the distance to a faraway star, etc. To measure such quantities, we find auxiliary easier-to-measure quantities x1,...,xn that are related to y by a known algorithm $y=f(x_1,...,x_n)$; then, we measure xi, and apply the algorithm f to the results $X_1,...,X_n$ of measuring xi. As a result, we get an estimate $Y=f(X_1,...,X_n)$ for $y$. This indirect measurement (data processing) is one of the main reasons why computers were invented in the first place, and one of the main uses of computers is scientific computing.

Measurements are never 100% accurate. The results Xi of direct measurements are, in general, different from the actual values xi. Therefore, the estimate $Y=f(X_1,...,X_n)$ is, in general, different from the actual (unknown) value $y=f(x_1,...,x_n)$. What do we know about the error $dy=Y-y$ of the indirect measurement?

In most cases, we know the upper bounds Di on the measurement errors $dx_i=X_i-x_i$ of direct measurements. Once we know such an upper bound, we can guarantee that the actual value xi lies in the interval $[X_i-D_i,X_i+D_i]$. In this case, the only information that we have about $y$ is that $y$ belongs to the range $f([X_1-D_1,X_1+D_1],...,[X_n-D_n,X_n+D_n])$.

Interval computations enable us to either compute this range exactly, or at least to provide an enclosure for this range. When the measurement errors are relatively small, we can expand f into Taylor series and ignore higher terms. If we only keep linear terms, we get an explicit formula for the range. For a quadratic approximation, all known methods of computing the exact range require $2^n$ steps - and since this problem is NP-hard, for large $n$, we can only compute enclosures.

In many cases, in addition to the upper bounds on dxi, we have partial information on the probabilities of different values of $dx_i$. In such cases, in addition to the interval range, we would like to compute the information about the probabilities of different values of $y$. There exist ways of extending interval techniques to such cases, see, e.g., S. Ferson, RAMAS Risk Calc 4.0, CRC Press, Boca Raton, FL, 2002.

One way to compute the enclosure of a quadratic approximation function f is to use naive (straightforward) interval computations. Due to the known dependence property, we often get excess width.

We do not have any excess width if the quadratic part is simply the sum of squares (i.e., if the corresponding matrix $A$ is diagonal). Every quadratic function can be represented in a similar form - as a linear combination of squares of eigenvectors. It is therefore reasonable to expect that (at least in some cases), this algebraic reformulation will decrease the excess width. Sometimes it decreases: e.g., if A has only one non-zero eigenvalue. In other cases - e.g., when all eigenvalues are equal - excess width increases.

In this talk, we analyze which method leads to a better estimate - by fixing eigenvalues and computing the expected excess width over all possible eigenvectors (with natural probability distribution). Result: eigenvector method is better $<->$ the variiance of eigenvalues is large.

Home page