Group Abstract Group Abstract

Message Boards Message Boards

2
|
1.1K Views
|
3 Replies
|
6 Total Likes
View groups...
Share
Share this post:

Problems with iterative process

Posted 4 months ago
POSTED BY: Denis Ivanov
3 Replies
Posted 4 months ago

So SetPrecision[tot, Infinity] helps a lot

POSTED BY: Denis Ivanov

Floating-point analysis

$\def\fl{\mathop{\text{fl}}}$

For roundoff error in double precision:

Let $\fl(u)$ represent the correctly rounded floating-point number closest to the real number $u$. Let $x,y,z$ be floating-point numbers between $0$ and $1$ (so that $x = \fl(x)$ etc.). Let $x+y+z$ be true sum and let $\fl(x+y+z)$ the computed sum. If $\fl(x,+y+z)=1$ (compared with or without tolerance), then $$x+y+z=1+\delta$$ where $\delta$ is some small roundoff error bounded by

$${-\tau_{-} \le \delta \le \tau_{+}}$$

with the tolerances $\tau_{-} \ge \varepsilon/4$, $\tau_{+} \ge \varepsilon/2$ and $\varepsilon = 2^{-52}$ (machine epsilon).

Now $2u$, is computed without roundoff error, as are $2u-1$ and $1-2u$, if $0.5 \le u < 1.0$. It is not hard to see that $x'+y'+z'$ will be either $1+2\delta$ or $1-2\delta$ (3rd or 4th rule resp.), if $(x',y',z')$ is the iterate following $(x,y,z)$. If $\delta \ne 0$, then eventually the iteration stops with "Total Error", unless a "Half Total" incidentally occurs earlier.

If $\delta=0$, then the following determines the behavior. If $0.5 \le u < 1$ and the trailing $n$ bits of the mantissa of $u$ are zero, then the trailing $n+1$ bits of the mantissa of $2u-1$ and $1-2u$ are zero. After at most 52 such transformations, $u = 0.5$ and you get a "Half Total". Since each of $x$, $y$, and $z$ would experience 52 such transformations if the iteration never stopped, a "Half Total" is inevitable.

Example: Case $\delta \ne 0$

seq = numsNest[{0.84474, 0.117199} // Append[#, 1 - Total[#]] &, 1.]
(*
{{0.84474, 0.117199, 0.038061}, {0.68948, 0.234398, 
  0.076122}, {0.37896, 0.468796, 0.152244}, {0.24208, 0.062408, 
  0.695512}, {0.48416, 0.124816, 0.391024}, {0.03168, 0.750368, 
  0.217952}, {0.06336, 0.500736, 0.435904}, {0.12672, 0.001472, 
  0.871808}, {0.25344, 0.002944, 
  0.743616}, {{0.25344, 0.002944, 0.743616}, "Error Total"}}
*)

The error in the total doubles at each step:

Total /@ Most[seq] - 1
Ratios[%]
(*
{0., -1.11022*10^-16, -2.22045*10^-16, 4.44089*10^-16, 
 8.88178*10^-16, -1.77636*10^-15, -3.55271*10^-15, -7.10543*10^-15, 
-1.42109*10^-14}

{ComplexInfinity, 2.`, -2.`, 2.`, -2.`, 2.`, 2.`, 2.`}
*)

The initial zero is because $\fl(x+y+z)=1$ exactly, but $(x+y+z)-1 = -2^{-54}$:

SetPrecision[{0.84474, 0.117199} // Append[#, 1 - Total[#]] &, Infinity] //
  Total // N[# - 1] &
(* -5.55112*10^-17  *)

The error after the first step is twice this, namely $-2^{53} = -1.11\times10^{-16}$.

Example: Case $\delta = 0$.

Remove the roundoff from an initial value less than 0.5 to get a total that equals 1:

SetPrecision[{0.84474, 0.117199 + 2^-54} // 
    Append[#, 1 - Total[#]] &, Infinity] //
  Total // N[# - 1] &
(*  0.  *)

seq2 = numsNest[{0.84474, 0.117199 + 2^-54} // Append[#, 1 - Total[#]] &, 1.]
(*
{{0.84474, 0.117199, 0.038061}, {0.68948, 0.234398, 0.076122}, 
  ... 
  {0.0136719, 0.244141, 0.742188}, {0.0273438, 0.488281, 0.484375}, 
  {0.945313, 0.0234375, 0.03125}, {0.890625, 0.046875, 0.0625}, 
  {0.78125, 0.09375, 0.125}, {0.5625, 0.1875, 0.25}, {0.125, 0.375, 0.5}, 
  {{0.125, 0.375, 0.5}, "Half Total"}}
*)

A visualization of the effective loss of precision:

MatrixPlot@RealDigits[#, 2][[All, 1]] & /@ 
  Transpose@seq2[[;; -2]] // GraphicsRow

enter image description here

POSTED BY: Michael Rogers
Posted 4 months ago

Thank you very much!
The problem is that by the results for rationals, there can be deterministic chaos.
And we need to understand that for reals can be natural behavior, and what - errors of rounding.
I see that the subject is not easy.

POSTED BY: Denis Ivanov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard