Re: Why is 1 smaller than 0?
- To: mathgroup at smc.vnet.net
- Subject: [mg73789] Re: [mg73749] Why is 1 smaller than 0?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 28 Feb 2007 04:38:46 -0500 (EST)
- References: <200702271049.FAA24070@smc.vnet.net> <45E455EB.1060305@wolfram.com> <20070228030819.GA10426@dirac.org>
Peter Jay Salzman wrote:
> On Tue 27 Feb 07, 10:01 AM, Daniel Lichtblau <danl at wolfram.com> said:
>
>>Peter Jay Salzman wrote:
>>
>>>This is an implementation of "steepest descent" to minimize a function. It
>>>runs for 36 iterations. On the 36th iteration, it claims:
>>>
>>> 0.0274 > .00001 True
>>>
>>>which is true. After the While loop exits, it claims:
>>>
>>> 0.03 > .00001 False
>>>
>>>which is true, not false, as Mathematica claims.
>>>
>>>
>>>Why is the While loop exiting prematurely? How do I write this so that it
>>>runs for as long as f(x,y) is greater than .00001?
>>>
>>>Thanks!
>>>Pete
>>>
>>>
>>>
>>>
>>>(* Implements Minimization via method of steepest descent. *)
>>>Clear[f, x, s, delf, a , xNew, iteration, delta, tolerance, theA, min];
>>>(* Function to minimize *)
>>>f[x_, y_] := 100`50*(y - x^2)^2 + (1 - x)^2;
>>>delf[x_, y_] := {-2*(1 - x) - 400.0`50*x(-x^2 + y), 200.0`50*(-x^2 + y)};
>>>
>>>
>>>iteration = 0.0`50;
>>>(* Initial guess. *)
>>>x = { {5.0`50,1.0`50} };
>>>(* Points "downhill" from the current position. *)
>>>s = {};
>>>
>>>While[ f[Last[x][[1]],Last[x][[2]]] > .00001,
>>>
>>> Print[f[Last[x][[1]],Last[x][[2]]], " ",
>>> f[Last[x][[1]],Last[x][[2]]] > .00001];
>>>
>>> (* Get direction to travel in (downhill) from grad f. *)
>>> s = Append[s, -delf[Last[x][[1]], Last[x][[2]]]];
>>>
>>> (* a tells us how far to travel. Need to minimize f to find it. *)
>>> xNew = Last[x] + a*delf[Last[x][[1]], Last[x][[2]]];
>>>
>>> (* Minimizes f with respect to a. *)
>>> {min, theA} = Minimize[f[xNew[[1]], xNew[[2]]], {a}];
>>>
>>> (* Update x using the direction *)
>>> xNew = xNew /. theA;
>>> delta = Norm[Last[x] - xNew];
>>> x = Append[x, xNew /. theA];
>>> ++iteration;
>>>];
>>>
>>>Print["Convergence in ", iteration, " iterations.\n",
>>> "Delta: ", delta, ", tolerance: ", N[tolerance], "\n",
>>> "Minimum point at ", Last[x], "\n",
>>> "Value of f at min point: ", f[Last[x][[1]], Last[x][[2]]], "\n"];
>>> Print[f[Last[x][[1]],Last[x][[2]]], " ",
>>> f[Last[x][[1]],Last[x][[2]]] > .00001];
>>
>>As someone pointed out, there is precision loss and at sufficiently low
>>precision values are treated as zero (hence less than .00001). To avoid
>>this pitfall you might work in any fixed precision that is no larger
>>than the input precision. This can be accomplished by wrapping, for example,
>>
>>Block[{$MinPrecision=40,$MaxPrecision=40},
>>...
>>]
>>
>>around your While loop.
>>
>>Daniel Lichtblau
>>Wolfram Research
>
>
>
> Hi Daniel,
>
> Thank you SO much for explaining how to get around the problem!
>
> I understand (loosely) that you're pinning the precision of the calculations
> and (even more loosely) why you use a Block[] instead of a Module[].
>
> I'd like to understand why I was losing precision in the first place. What
> are the rules that govern precision loss during a calculation?
>
> (You sent this via private email, and it's bad form to forward private email
> to a list, but if you think this would be a good thing to put into the
> Mathgroup archives, feel free to respond via the list. Seems like a good
> thing to post publically).
>
> Many thanks!
> Pete
Okay, I'll respond to MathGroup and make a few comments on questions and
possible puzzlement.
(1) Yes, I'm pinnint the precision to some set value. I think there are
ways to mess this up e.g. by starting with machine numbers, but that is
not an issue in your example.
(2) Block will indeed block the effect of the global
$MinPrecision/$MaxPrecision. Why do you think they call it "Block"?
Well, okay, the real reason for the name is that it is a scoping
construct, remotely similar to the notion of a code block in more
classical languages. But it just happens that the name fits this
variable scoping/masking quite well. Module (which of course is also a
"block-like" code construct) instead creates its own local versions of
variables, hence will not affect globally understood system variables as
we wish to do here.
(3) Precision loss is part of the game with significance or interval
arithmetic. Your iteration will use numbers beginning at 48 or so
digits precision and in the process of doing arithmetic upon them (say,
to compute f[] or delf[]), first order approximations of error are made.
These in turn contribute to the "fuzzing", that is, loss of precision.
Also you may be losing it in other ways e.g. in "black box" functions
such as Minimize. But I don't think that'a an issue in your specific
code. Actually I thik most of it happens in delf.
(4) Another way to work around the problem is to use SetPrecision to
arbitrarily increase precision of, say, the delf result (and maybe other
computed quantities) in each iteration. Since you are doing an iterative
convergent process, this will not cause trouble provided you do not try
to recover a result with smaller residual than the original input
accuracy/precision will allow. With your epsilon of 10^(-5) this should
not pose a problem.
Daniel Lichtblau
Wolfram Research
- References:
- Why is 1 smaller than 0?
- From: p@dirac.org (Peter Jay Salzman)
- Why is 1 smaller than 0?