Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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




  • Prev by Date: Re: Re: conditional is giving wrong value
  • Previous by thread: Why is 1 smaller than 0?
  • Next by thread: Re: Why is 1 smaller than 0?