Accuracy and error analysis (was Re: [Maxima] primes)

[C Y , Sat, 14 May 2005 09:33:26 -0700 (PDT)]:
> True.  From the speed point of view though, I'm less worried since
> there are different levels to choose from.

Well, typically it goes like that: You work with some theory,
implement the relations you need in the CAS in order to have something
to play around.  Then you want to compare your results with some data
sets, and you say to yourself: I already have this implementation of
the theory, let's just use that.  When the data sets are large, you
then easily run into speed problems: It is then really important to
get things to run as fast as you can manage, in particular if you want
to maintain the interactivity of a REPL rather than to batch
everything.

For trivial things there are never problems, I guess, and things are
fast enough most of the time anyway.  But the converse case is
sufficiently frequent in my experience.

(Whether it is a good idea to analyze large sets of numerical data in
a CAS is a different question altogether.  I would guess that Maxima
is in a better position in this respect as you can always fall back to
the underlying Lisp if need be.)

> if you care about significance you HAVE to pay those costs, if you
> want a reliable answer.

Caring for the significance of the answer is something different from
trying to evaluate arbitrary expressions to a given precision.

In current Maxima, things are pretty clear: There is no uncertainty
attached to the numbers, numerical computations simply give you
whatever the underlying Lisp happens to give you.  Things are simple,
but the system does not give you any help in figuring out how much of
this is to be trusted, and the user does not have any reason for
trusting the CAS to tell her something about the uncertainties.

One may also do some (semi-automatic?) error analysis: These
quantities are known with those uncertainties, and the result of some
evaluation then also has some uncertainty.  Provided things are
implemented correctly, this analysis can be trusted, and everything is
well defined.

But the question "what are the first N digits of f(x), where x is only
known to M significant digits" is, I think, not well posed in general.
Certainly what Mma does in this respect is quite problematic, v.i..

OTOH, the answer to the related question "how large may the
uncertainty in x be so that f(x) has N significant digits." would
again be well defined (and intimately related to the automated error
analysis) but, I am afraid, less appealing to most people.

> I'm afraid I don't understand the necessity of extra care in
> defining functions

I don't know how exactly you want to handle things, so let me clarify
the Mma way which is what I was referring to:

Suppose you are computing f(g(x)), where x is known to within dx.  You
want to get a certain accuracy for the result.  So you can start with
x+-dx, compute g(x+-dx), and feed that result to f.  It may now happen
that you do not get a result of the desired accuracy.  What do you do?

Mma's strategy is to reduce dx and repeat the whole calculation.  And,
in fact, to do so as often as is needed to fulfill the promise "I am
going to give you as many significant digits as you asked for" (which
often masks the loss of significance in its significance arithmetic,
at a cost in the number of digits carried around).  Only once some
limit ($MaxExtraPrecision, if memory serves right) is reached, does it
give up and give you a lower-precision result.

This way of doing things has at least three problems:

- if dx was meaningful from the outset, re-setting it to something
  smaller is not valid.

- it is not clear when to give up: if f is lambda([x],
  constant-of-fixed-accuracy), you cannot raise the accuracy of the
  result.

- the number of evaluations of f or g can no longer be predicted.
  Which is fine when those functions are pure as they should be, but
  this is not not necessarily so in the presence of side effects.  And
  unfortunately people used to imperative languages like to write that
  kind of code even in Mma.

  Is this a problem in Maxima?  I don't know, it depends on the
  population of users it attracts.  In Mma it is a problem, as
  attested to by all the bugs that come from users not understanding
  Block[] and Module[] and With[] etc. and wildly use
  CompoundStatement[] instead; in Maxima, that corresponds to
  forgetting about, e.g., block() and using (expr1, expr2, ...)  all
  the time.


  In current Maxima I am not aware of situations where it is difficult
  to predict how often a function will run.  Maybe I am missing
  something - there certainly are some things I have never or hardly
  ever used.  This means, in particular, that side effects in
  functions are not necessarily bad even when the side effects affect
  the value the function returns.
  
  In particular,  in the current scheme the ugly code
  
      foo(x) := (myPrivateVariable: myPrivateVariable+1,
                 bar(myPrivateVariable,x)) $
  
      block([myPrivateVariable: 0],
            thru 10 do something(foo(1.)));
  
  will work.  If, however, give_me_to_N_digits() follows a strategy
  where computations are repeated with higher precision,
  
      block([myPrivateVariable: 0],
            thru 10 do something(give_me_to_N_digits(foo(1.))));
  
  will be a nasty surprise.

Regards,

Albert.