Stability in scientific computing

Goals:
- Provide approximations to mathematical formulations
- e.g. numerical integration, root finding, matrix computation,
differential equations
- Accuracy: how close the approximation is to the truth
- Efficiency: how quickly the answer is obtained
- Stability: how "well-behaved" the approximation is

Historical perspective:
- Newton (calculus was only approximate)
- Gauss (linear systems)
- lots of others
- Turing: war effort, build a digital computer (w/ Wilkinson)
- von Neumann: US computer -- apply to particular scientific problems
(especially codes, turbulence)

Fundamental issues:
- Floating point computation (ring is no longer associative,
arithmetic not exact)
  - Can the small errors in rounding lead to arbitrarily bad
  approximations? 
  - Early estimates (Hotelling, von Neumann, Turing) for Gaussian
  elimination implied up possible growth in roundoff error
  - Some methods become notoriously unstable.
  - Need stability of methods under small perturbations
- Devising mathematically stable algorithms
  - Look at model partial differential equations and what can
  go wrong with simple finite difference methods
  - von Neumann stability analysis
- Overcoming algorithmic implications of ill-conditioning
  - Poisson equation in 1d: condition number
  - Summary of various algorithms:
    - Dense GE
    - Sparse GE (but what about 2d, 3d?)
    - Jacobi/Gauss-Seidel/SOR
    - Krylov
    - MG


Floating point computations:
- Example: evaluating a polynomial:
  - Naive way
  - Horner's method
  - What if you have the roots?
- What are implications for, say, finding the root by bisection?
  - Could terminate anywhere in a relatively large interval!
- Floating point model:
  - sign plus exponent plus fraction (1/8/23 single vs. 1/11/52 double)
  - basic laws of operation
  - cancellation: roots of quadratic equation
    - b^2 - 4ac, what if b^2 >> |4ac|, cancellation (take sign(b) for
      larger root)
    - b^2 \approx |4ac| --> must use extended precision, or rewrite problem
       - b^2 = 100, 4ac = 99.99 (= 100 in 3-digit arithmetic)
       - For any precision, can break it (but must be more perverse)
  - forward error: what is error in | y - y* |, where y = f(x), y* = alg(x)
- Backward error (statement about the algorithm, not the problem)
  - Backward analysis for Horner
    p = a_d
    for i=d-1,0
      p = x p + a_i

    call each stage of computation p_i so we can do analysis on it

    p_d = a_d
    for i=d-1,0
      p_i = x p_{i+1} + a_i
    p is p_0

    p_d = a_d
    for i=d-1,0
      p_i = ( (x p_{i+1})(1+\delta_{1i}) + a_i)(1+\delta_{2i})
    
    p_0 = \sum_{i=0}^{d-1}
    [ (1+\delta_{1i})\product_{j=0}^{i-1}(1+\delta_{1j})(1+\delta_{2j})]
      (a_i x^i) + 
    [ \prod_{j=0}^{d-1}(1+\delta_{1j})(1+\delta_{2j}) ] a_d x^d 

    don't keep track of every detail, note that:
    1-j\eps \leq (1+\delta_{i1})...(1+\delta_{ij}) \leq (1+\eps)^j \leq 1+j\eps

    as long as j\eps << 1.

    p_0 = \sum_{i=1}^d (1+\bar{\delta}_i)a_i x^i, |\bar{\delta}| \leq 2d\eps

    = \sum_{i=0}^{d}\bar{a}_i x^i

  - roundoff error causes you to solve (exactly) a nearby problem in
    the sense of perturbed input data (Wilkinson's Turing award)
  - couples to perturbation theory for the problem.  How sensitive is
    the problem to small changes (e.g. roundoff)?

- Conditioning
  - Relative condition number of a function:
    	    | f(x + dx) - f(x) | / | f(x) |
	      	    \approx
	    | dx | / | x | * | f'(x) |  | x | / | f(x) |
	    relative: | f'(x) | | x | / | f(x) |
	    absolute: | f'(x) |
  - Backward stability:
    error = | alg(x) - f(x) | = | f(x+dx) - f(x) | \approx |f'(x)| |dx|
    --> only "ill-conditioned" problems can have large absolute error!
  - Condition number for polynomial by Horner's method

    | p_0 - p(x) | = | \sum_{i=0}^d \bar{\delta}_ix^i |
       \leq 2d\eps \sum_{i=0}^d  |a_i x^i|

    If roundoff error goes \bar{delta}_i = \eps \sign(a_i x^i), we
    get within 2d times the error bound

    Relative condition number of polynomial p = \sum a_i x^i
    k(p) \sum_{i=0}^d |a_i x^i| / | \sum_{i=0}^d a_i x^i |

    High relative error in computing things close to p(x) = 0,
    infinite relative error in computing p(x) = 0.

  - Condition number relates to "distance to singularity"
  
    p(z),q(z) polys with coeffs a,b:

    define d(p,q) = \min_\alpha : |a_i - b_i| \leq \alpha |a_i|
    if a_i \neq 0, then d(p,q) = \max_i |a_i - b_i| / |a_i|

    Theorem : p(z) is not the constant 0
    then min(d(p,q) : q(x) = 0) is 1 / k(p)

  - Suppose you know the roots...Beats condition number problem
    p = p(x)(1+\delta) for |\delta| \leq 2d\eps

- Linear Algebra
  - What's the norm of a matrix? \sup_{||x||=1} || A x ||
  - Condition number of a matrix ( k(A) = || A || || inv(A) || )
  - Backward analysis for Gaussian elimination (tedious)
    - only get it for partial or complete pivoting
    - analysis is wildly unstable without pivoting!
  - Perturbation theory for linear systems says 
    Ax = b ;; (A + dA)(x + dx) = b + db
    dx = A^{-1} (-dA xhat + db)
    || dx || / || x ||
      \leq 
    k(A) ( || dA || / || A || + || db || / ( ||A|| || xhat|| ) )    

Bibliography:
http://docs.sun.com/source/806-3568/ncg_goldberg.html (What every
computer scientist should know about floating-point arithmetic)

Demmel, Applied Numerical Linear Algebra, SIAM 1997.

Higham, Accuracy and Stability of Numerical Algorithms, SIAM 2002.

Wilkinson, Rounding Errors in Algeraic Processes, Dover 1994.