Re: Why Mathematica does not issue a warning when the calculations
- To: mathgroup at smc.vnet.net
- Subject: [mg117795] Re: Why Mathematica does not issue a warning when the calculations
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 31 Mar 2011 04:04:52 -0500 (EST)
Richard Fateman wrote:
> On 3/29/2011 3:58 PM, Daniel Lichtblau wrote:
>> .. snip..
>
>>
>> (DL) 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.
> (RJF) Considering this other system (sorry, Steve C, will you bear with
> us?)
>
> Let z:= -251942729309018542,
> a:=evalf[10])(z),
> b:=evalf[20](z),
> for val10, you computed f(a)
> and for val20 you computed f(b),
> and certainly the results are different.
>
> But contrary to your contention, this particular system does NOT
> believe that a and b are equal. In particular,
> a-b is .9018542 x10^7, and the test
>
> evalb(a=b) returns false.
>
> Your test
>
> type(a=b,equation) returns true.
>
> But this does not test for equality. It is a kind of pattern match that
> confirms the type, namely an equation.
> In particular, type(3=4,equation) also returns true, because the TYPE
> is "equation".
>
> Translated into Mathematica, this is approximately MatchQ[foo == bar ,
> HoldPattern[Equal[__]]]
Good catch.
I hereby confirm that my example does not apply to that unstated other
program (band plays: "Smoking in the Boys' Room"), insofar as a!=b in
that system, for that problem. I will later provide a more contorted
variant where the pathology does appear.
I will skip some of what you later wrote, but state here that you are
correct: the truncation issue arose earlier than I had claimed.
Technical remark: It still is "truncation error", not in the sense of
being wrong (it isn't) but in the sense of showing why it is not in
agreement with an exact computation. This does not argue for or against
either one of us.
>> 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.
> Here is a fundamental disagreement. When a human types 0.1 as a number
> to be represented as a hardware float (binary) it is converted to binary
> and has a certain bit representation, B. That now is the number that
> the computer has in its electronic memory. B is close to, but not
> exactly 0.1 decimal. It is closer to 0.10000000000000000555 except
> rounded to "machine precision". it is 0.1.
> People of a certain sophistication will realize that not all decimal
> numbers have exact binary representations. Just the fraction 1/3 does
> not have an exact decimal representation. We must live with some
> version of this radix stuff and hope we can explain it when necessary.
Yes, but...I will surmise that any program using binary representation
and decimal input will exhibit the pathology my example purported to show.
For a specific example, I can get the bad behavior by round-trip
conversion to-from binary.
> Digits := 30:
> f := x -> evalf(convert(x,decimal,2) - (6 + sqrt(2))^20/10^6):
> v15 := convert(evalf[15](251942729309.018542),binary,60):
> v30 := convert(evalf[30](251942729309.018542),binary,60):
> evalb(v15=v30);
true
> val15 := f(v15);
val15 := 0.000457127208979235
> val30 := f(v30);
-6
val30 := -0.874881474866 10
> evalb(val2=val3);
false
A bit contorted, I admit. Mayber outside the scope of general usage. But
it serves to show specifically that starting with decimal and having
binary under the hood will not, in general, play nice with silent
appending of zeros, whether digits or bits, past the last significant
figure to the right of the radix.
> According to Mathematica, B, in binary, can be written out as something
> like
> 0.0001100110011001100110011001100110011001100110011001101... base 2.
> If we add binary zeros to the right of B, we do not change the value. B
> is not muddied by adding zeros. The number in the computer is just as
> accurate, exactly as accurate, as B, as an approximation to 0.1, even
> with more zeros.
Yes, but... if you have truncated in one case (hence losing those later
bits) but not the other, adding zeroes will pollute the one and not the
other. See above example. (Which I really, really hope is correct this
time around. You cannot imagine the scorn I will face if I lose the same
debate TWICE to a retired professor.)
> I, on the other hand, see no reason to assume that some number in the
> computer, is not a perfectly valid number, exactly what the user wants
> to compute with, and is accurate to all digits and can be extended by
> zeros if need be.
> [snip...]
> Executive summary of the paragraph above:
> You couldn't possibly mean THAT number, which is so close to 1/10. So
> I will refuse to allow you to do arithmetic with it.
Nothing wrong with what you wrote, just a different view from my own.
I'll give an example in a separate post.
> .. snip..
>> 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...).
> Can you cite a reference for this? This is the essence of significance
> arithmetic, but really not prevalent.
It's fairly necessary in bigfloat linear algebra. I believe it is also
needed in other "standard" numeric algorithms/functions such as
quadrature, root finding and numeric optimization.
But yeah, this is all belief on my part and I do not have references. My
upcoming example might give an idea of why I hold to such (apparently
unusual) ideas.
>> This is because it means we do not pollute the arithmetic with
>> manufactured digits (which might arise if we do differently and raise
>> precisions).
> Most other arithmetic systems attempt to retain as much information as
> possible. You claim that extra precision is "manufactured digits",
> though you (perhaps) allow extra digits that happen to correspond to
> decimal digits.. e.g.
> 0.25 decimal is 0.01 base 2. exactly. So 0.010000000000000 base 2 is
> exactly 0.25 decimal. Are those digits better?
> I claim that these digits (binary zero digits) can be useful. Sometimes
> they allow numbers to cancel WHEN THEY SHOULD, and leaving them out
> causes more problems.
They can be useful when there is no prior truncation error in the
construction of the number. I don't think I ever stated or implied
otherwise. But again, extension by zeroes implies knowledge that it is a
vlid thing to do (so don't do it blindly). Again, consider the case
where we start with .3 rather than .25. If you think, "Well, obviously,
you can rationalize and then numericize to higher precision", then
conside instead .3217784, and tell me why I should believe
rationalizing/numericizing is a good thing to do for that.
> Nevertheless, a reasonable design for a system is where one judges the
> precision of each operation by the precision of the target. Presumably
> it has to do some work to determine how much precision is needed of the
> operands -- extended if necessary by zeros, to achieve the desired
> target precision.
> ... snip...
That works fine for many types of numeric computation. Not all, though.
And for some problems, the conditioning of the problem is such that
attaining the precision of, say, the most precise operand, is simply not
possible. For others, e.g. root-finding, I don't think it is even
well-posed how to do that. For example, how do you know your root is
precise to x digits, if your input is only precise to x/2 digits?
Sometimes you can check residuals, do Newton iterations, etc. But
eventually your steps only produce fuzz because input is not
sufficiently precise to get meaningful digits.
> Here's an approximation to Pi. 3.1 decimal
> There is absolutely nothing wrong with taking that number (31/10) and
> writing it out, in decimal or binary to any number of places.
> It does not become more accurate as an approximation for Pi, just more
> precise as an approximation for 31/10.
Yes. But doing that will not give a more accurate result in a
computation, if the original value was in fact Pi.
>> 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
>
> OK, this, above, is the number that the computer knows about.
>>
>> 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.
> Sure. That's the closest binary number to the one you typed in.
Well, not really. Here you are tacitly (I think) inferring precision of
input. One can get arbitrarily close in binary by going to a higher
precision representation.
What I mean:
In[4]:= approx13bigprec = FromDigits[RealDigits[1.3`25,2],2]
Out[4]//InputForm= 6286414261996071708472115/4835703278458516698824704
In[6]:= N[approx13bigprec,25]
Out[6]//InputForm= 1.29999999999999999999999995864096937235`25.
>>>> 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.
I agree with that. But one does not always work with integers. Or with
"nice" decimals (defined howsoever you like, provided 1.2345678765998821
is not regarded as "nice").
>> 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.
> I think the dearth of papers on significance arithmetic is some
> indicator of popularity.
Maybe. An alternate view, common in science, is this. One writes about
the things for any or all of these reasons: one understands them, can
analyze them, can create/implement/work with them in a lab
(computational, in our case) setting, has prior experience handling
them, can get the article past peer review, etc. This might or might not
be applicable in the case under consideration.
>> Ah. Though I wonder if anyone else is still awake?
>>
> My fans? Your fans? people who use Google groups search for error
> analysis, significance arithmetic, your name, my name...
My name might attract attention. But only because there are (at least)
three of us, and one is in Hollywood.
>> 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.
> um, not really.
Well, actually yes, at least until you find a flaw in the upgraded
example. Which possibility I do not rule out, given the sophomoronic
mistake of its predecessor. (Yes, that's really a word. For today.)
>> 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 would require not only identical hardware and computation, but
> identical status settings for traps etc. Also, no hardware errors (such
> as might be caused by cosmic rays, radioactive decay of elements in the
> computer chips, mechanical breakage etc.).
> I am not aware of a circumstance in which different results occur.
> Sounds like one of those riddles.
Here are two ways it might, and in practice does, happen.
(1) In some machines, double precision array alignment affects whether
one type of register is used or another. In practice I believe this is
double or extended, that is, 64 or 80 bit. Different results can and do
emerge. Numerical analyts are fine with this, and so am I. Just showing
that it can and will happen in surprising settings such as identical
IEEE-conformant computations on the same machine.
(2) Related, but different, and probably not conforming to your added
requirements: Optimization level of a compiler will affect whether or
when a register is stored to memory. On many machines double precision
computations are performed in extended precision registers. The store
will truncate. In addition to optimization level, swapping out of a task
in a multiple processor environment can also make this come about.
(2b) I don't think that truncation is a bad thing, but it already gives
rise to the phenomenon under discussion (same hardware, different
results). But far worse, extended prec registers have both larger
mantissae AND exponents. Rationale: if there was an intermediate
overflow then having larger exponents will often allow the computation
to move forward). So here is the problem: a a store to memory at a bad
time will cause an overflow.
For item (2) I cannot tell you how annoying it is to track down such a
problem, as a debug version of a program will behave differently from a
production one due to different compiler optimization settings.
There is a nice discussion of item (1).
http://thenumericalalgorithmsgroup.blogspot.com/2011/02/wandering-precision.html
>> 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.
> Judiciously choosing among conflicting requirements is sometimes an
> important element of design.
>
> RJF
I fully agree.
Clearly you and I have different preferred resolutions. But that is
mostly a matter of having differing priorities as well as coming into
this with different backgrounds and expectations. Simply stated, we are
quibbling over minor resolution details, such as whose numerics mojo is
more powerful.
(one of three) Daniel Lichtblau
Wolfram Research