Re: Why Mathematica does not issue a warning when the calculations
- To: mathgroup at smc.vnet.net
- Subject: [mg117725] Re: Why Mathematica does not issue a warning when the calculations
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 30 Mar 2011 04:14:24 -0500 (EST)
> But then it is clear that you are testing a quality of x that varies
> with the precision of
> x, so that even if x1==x2, their precisions may differ, and so the
> results may differ.
This was exactly what you had requested, specifically, a case where this
would happen using fixed precision arithmetic.
> I think the issue you are trying to raise is that any system with
> variable precision must decide on some precision to use to evaluate
> constants like sqrt(2). If that precision is derived from the precision
> of a nearby float (not its value), then that constant, or the expression
> involving it, may come out differently for different precisions. That's
> not a requirement, though.
Here is an example I trust you will recognize.
f := x -> x + (evalf[40](6) + sqrt(evalf[40](2)))^20:
Digits := 40:
> val10 := f(evalf[10](-251942729309018542));
val10 := 9018542.6877782033820121476772
> val20 := f(evalf[20](-251942729309018542));
val20 := .6877782033820121476772
You will notice that the constants are all evaluated numerically in
advance, and all to 40 digits.
> Aha, I finally see your point here. You think that by using With
> [{prec=10 }....] you are doing what non-Mathematica programs do,
> silently. That is, you believe they scan all the numbers in the
> expression and find the precision with which to evaluate the Sqrt[2].
> Whether to use the lowest or highest precision, who knows. In this case
> there is only one precision and so it doesn't matter. What is essential
> here, and which I now see I missed, ( because it is inconsistent with
> how some other systems actually work), is your assumption that the
> precision of the input would be determined in this way.
The example I show above indicates that one can set the precision of
that square root quite high, and still get a "bad" result. This is
truncation error (of the argument to f) followed by cancellation error.
How high precision is used for the other stuff does not really matter.
In an example of this sort a program COULD, in principle, raise the
precision of the input argument. That will only work if that input is
exactly representable in the underlying arithmetic. For decimals, in
programs that use base 10 behind the scenes, this can be done. If the
input is decimal and the representation is binary and the decimal is not
exactly representable (e.g. 1.3 in binary) then you'd be in trouble.
Likewise if the input was a numericized approximation to pi, say.
>>> OK, you provided an example that shows a == b, f[a]=!=f[b], in
>>> Mathematica.
>> I provided considerably more.
>>
>> I showed that this happens with fixed precision. So you will
>> see it with other programs as well, provided they support an
>> equality test that allows for
>> lower precision bignum == higher precision bignum
>> I moreover showed that at least one such program has that
>> behavior.
>>
> You showed this for the case where the evaluation of an expression is
> done differently for different precisions because the evaluation
> mechanism extracts from the inputs, an extra parameter, namely
> "precision" and does a different calculation.
Again, this is not really essential to the process in general (as
indicated in example near top). Though it is necessary if we want
Mathematica to emulate fixed precisiona rithmetic in any sensible way. I
do not regard manufacture of digits or bits, that is, raising precision
artificially, as "sensible", unless one has a priori a justification for
doing so.
> If you had said so in
> the first place it would have been clearer (to me). If the system's
> evaluation program determines the desired precision by some other
> mechanism, something that is actually pretty likely, then your
> simulation in Mathematica is simulation of something that isn't
> pertinent.
I believe the implicit claim that, in principle, the assigning of
precision to exact values might make a difference.
In practice, and not just in Mathematica, a coercion model that forces
arguments to the lowest precision in the expression is a sensible way to
do computations (there are exceptions...). This is because it means we
do not pollute the arithmetic with manufactured digits (which might
arise if we do differently and raise precisions).
> For example is there a global "precision" (this is
> common) that governs the calculations? Or is the precision is
> determined by the target location as for an assignment z=f[x], look at
> the precision of the "container" for z. This is the method used by
> MPFR, which I think Mathematica uses internally.
Mathematica uses GMP for integer arithmetic. MPIR is a fork of GMP, not
used in Mathematica. Neither does Mathematica use MPFR for bigfloat
arithmetic.
Any or all of this could change in the future. But to date neither MPIR
nor MPFR have been used under the hood in Mathematica.
>>> (though the two functions are different...) But your effort agrees
>>> with
>>> my point, that a==b does not mean f[a]==f[b] in Mathematica.
>> Or several other programs, I dare say. Did you try any of them
>> on the example I showed?
>>
> Um, try them using some other computer algebra system? For those
> numbers, a is not equal to b, so the example doesn't work.
Actually I showed quite explicitly that they are equal in one such
program. To repeat:
> type(evalf[10](-251942729309018542) = \
> evalf[20](-251942729309018542), 'equation');
true
>> Not so. I contend that any program that claims 2.000==2.0000000
>> will do what I showed.
> You, as well as I ,fall into thinking inside a box. You assume that the
> precision of the computation of an expression is derived from the
> (lowest?) precision of any of the constituent parts. I assume it is not.
Okay, this could make a difference. The lines below shows that at least
one other program must do something like a coercion to lowest precision.
> type(evalf[10](21/10) = evalf[20](21/10), 'equation');
true
(And more compelling...)
> type(evalf[10](21/10) = evalf[20](2100\
> 0000000001/10000000000001), 'equation');
true
I do agree that if coercion to lowest is not done, then outcomes like
those above are unlikely. I also happen to think that would be a bad
outcome.
>>> OK, here is what I think SHOULD happen. I use a hypothetical computer
>>> algebra system that resembles Mathematica but not entirely...
>>> [snip...]
>>> a-b --> convert both to longest precision and subtract.
>> Why does this make sense? You are manufacturing digits (which
>> happen to be zero) for the lower precision value.
> Yep.
>> This will
>> do Really Bad Things if the underlying representation is
>> binary-based.
> Nope. You add binary zeros.
To spell out my thought in more detail, it will do bad things if the
input was decimal and not exactly representable in binary, e.g. 1.3.
Likewise if the input was, say, pi and had been numerically
approximated. This would (often) be bad because extending with zeroes
would not recover the value in question.
Let me show what I mean in regard to 1.3. If we take the bit pattern and
recover the "exact" value, we
In[10]:= approx13 = FromDigits[RealDigits[1.3, 2], 2]
Out[10]= 5854679515581645/4503599627370496
In[11]:= N[approx13, 22]
Out[11]= 1.300000000000000044409
This is what would be the result, were we to add zero bits to the
(truncated) binary representation of 1.3.
The only way to avoid this would be to work entirely in binary,
including input/output. That way you only have exactly representable
values. But again, numerical approximations to things like pi would
still not behave very well.
>> In effect you'd be putting garbage right where
>> you are trying to avoid it.
> Nope. Consider representing the integer 10 as a float. Do you mind
> adding decimal zeros or binary zeros to get 10.0000000000000000 ?
> Are you adding garbage? I think not.
And I agree, for this type of example. But not all numbers will have
exact representations in the underlying base. For those, putting zeros
at the end to raise precision would just (as in example above) be
putting junk into the number beyond its original precision.
> Say that you have the
> single-float 10.0. Do you think that it represents the same number as
> the integer 10? Or do you think it represents 10.000000xyzw .. where
> xyzw could be any digits at all? Now if that 10.0 were, say, a
> measurement that you knew was approximate, ok. But a computationally
> more useful model (in the opinion of most people) is that the
> single-float 10.0 is a compact representation of exactly the real number
> 10, that happens to have an exact representation in decimal or binary
> floating-point format.
> This model allows you to do careful analysis of errors etc. You can
> even simulate the other model with xyzw fuzz, using this.
Actually I take this approach for computing numeric Groebner bases. But
the rationale is quite different. The pretense, as it were, is that I am
solving a nearby problem. For most purposes, where such bases are used,
that is fine.
I'd be a bit surprised if your model is really the preferred one in
either the numerical analysis or scientific computation communities. But
I've not done any sort of poll.
> It is much harder to simulate the exact representation idea if the
> arithmetic you are presented with has the xyzw fuzz already, and you
> have to do gyrations to "get the fuzz out".
Yes, maybe, I suppose...but the fuzz (or badness) is already there if
your input was not exactly representable. This is true for your model as
well as any other.
>> I'm not sure that adding binary (as opposed to decimal) zeros
>> would be any better an idea, in this setting.
> It is the same for integers. For fractions it is also about the same.
> Unless you are wedded to the necessity of exact representation of
> certain decimal quantities such as 0.01. This has an advantage for
> bankers, and especially in inflated currencies or long-running mortgage
> payment calculations.
I'm not. i simply point out that adding bits to 0.01 will only make it
look more precise. It won't be (in the standard nomenclature) more accurate.
>>> now define the function f[x,newprec]:= N[x+(6+sqrt(2))^20 ,newprec]
>>>
>>> f[b,n] == f[exact,n] FOR ALL n, because b==exact in this system.
>> This will in general not be achievable, if n is larger than the
>> precision of b. At least not with semantics that does comparison
>> by raising rather than lowering precision.
> In this case b and exact represent the same integer, b in a
> floating-point format, exact is in an integer format. This should
> assure that f computes the same value.
Seems workable, if b is in some sense known exactly (e.g. is a numeric
approximation of an integer). This sort of input is seen in SOME real
problems. it is by no means what we have in ALL real problems.
>> (DL)
>>>> I think what you want is a different language...
>> (RJF)
>>> Or the same language with different semantics.
>> Which would be sufficiently alien to Mathematica as to make it a
>> different language.
> Actually, I doubt it. Indeed, if one were able to drop into Mathematica
> a different arithmetic model, how much of the test suite would run
> exactly the same? Of course those parts that explicitly rely on floats
> would mostly not do the same thing. And there may be parts that
> implicitly use floats such as heuristics to see if something is zero for
> simplification purposes. But my guess is that lots of stuff would work
> the same.
Numeric zero testing would go belly up. In consequence...it would take
with it some amount of other things in the exact computation realm.
Limit and Series, for example, would be in trouble (hence also definite
integration). I believe some aspects of Reduce would be slowed beyond
usability. Exact linear algebra would suffer, though more modestly.
In the numeric arena, NSolve would be dead.
Functions that rely on iterative processes, such as FindRoot, would most
likely be fine. Likewise most (though not all) numeric transcendental
function evaluations would work. Exceptions might be those that need to
"round trip" to determine appropriate precision. Although even those
might still usually work just fine.
That said, you are probably correct that most tests would pass (well,
not the ones that are already failing...). But the changes to those
items noted above would make for what I'd regard as a very different
Mathematica.
> If someone raises a question about computer algebra systems in the
> context of Mathematica, and there seems to be a possibility that I can
> contribute to educating the poster and others reading the posting about
> that issue, then I sometimes chime in. Sometimes the questions are
> "why does Mathematica do X" and the answers from others look like
> "Because it does X" or "It has to do X" or "Everyone knows that X"
> or "you should ask Y". But my reaction might be more along the lines
> of "Mathematica was designed/programmed to do X, though it is not
> inevitable that a computer system should do this. I can understand why
> you might expect, from mathematical knowledge and intuition, that
> Mathematica should do something else, and it would be nice if it
> fulfilled your expectations, and here's how it might be changed to do
> that."
> Or something along those lines.
Ah. Though I wonder if anyone else is still awake?
>> I'd be more charitably disposed, had you simply recognized from the
>> beginning that my example met your request.
> Ah, but perhaps you now understand my point -- you thought that your
> programs, using With [ {prec..}..] correctly modeled the behavior of
> other programs. I, on the other hand assumed quite the opposite, that
> you were inserting two different constants in a program, and then you
> contended that they were the same! I probably had as much concern for
> your sanity as you for mine. How could two textually different programs
> with different numbers in them be the same? I do see (now) what you
> were driving at.
I leave open the issue of how a variant of my fixed precision code,
properly ported to an arbitrary different program, might or might not
conform to your requested counterexample. Clearly at least one such
program behaves as I expected. I rather expect a program withn the
semantics you have in mind would fail to provide a counterexample by
flunking the a==b criterion. Well, depends on how equality is tested,
(slightly independent of) whether precision is raised or lowered.
The mention of bitwise 'sameness", as opposed to some program's idea of
"equal", jarred my memory. Maybe you know this, but for the sake of the
remaining third reader (if still breathing) I will state it. One can have:
(1) Identical machine double input (down to the last bit)
(2) Identical harware
(3) An identical arithmetic computation
yet still obtain
(3) Different results.
What's more, this is all quite within the realm of IEEE compliance.
I don't know whether or how this fits into the numerics you would like
to see. Just wanted to point out that consistency might be deceptively
hard to attain.
Daniel Lichtblau
Wolfram Research