Message Boards

WOLFRAM COMMUNITY

115 Views

3 Replies

2 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematics Mathematica Calculus Numerical Computation

Problems with iterative process

Denis Ivanov

Posted 1 day ago

For some reasons I need to consider the following iterative process. We have three positive numbers with constant total such that $x+y+z == t$ Consider the iterative algorithm: if the sum of numbers is not equal to $t$ - stop (with message "Total error") if one of the numbers is $t/2$ - stop (with message "Half") if one of the numbers is greater than $t/2$ (let it be $z$) - returns new numbers $2z-t,2x,2y$ if all numbers are less than $t/2$ - return the new numbers $t-2z,t-2x,t-2y$ ClearAll[genNumsReals]; genNumsReals[ tot_] := Block[{a = RandomReal[{0., tot}], b}, b = RandomReal[{0., tot - a}]; {a, b, tot - a - b}]; ClearAll[genNumsRationals]; genNumsRationals[e_, tot_] := Block[{a = Rationalize[RandomReal[{e, tot - e}], e], b}, b = Rationalize[RandomReal[{e, tot - a - e}], e]; {a, b, tot - a - b}]; Clear[nextNums]; nextNums[nums_, tot_] := Which[ Unequal[Total@nums, tot], {nums, "Error Total"}, MemberQ[nums, tot/2], {nums, "Half Total"}, AnyTrue[nums, GreaterThan[tot/2]], Mod[2 nums, tot], True, tot - 2 nums]; numsNest[init_, tot_, iterMax_ : 100] := NestWhileList[nextNums[#, tot] &, init, (UnsameQ[##] && Length@Last[{##}] == 3) &, All, iterMax]; Firstly, is it equal to `Mod[-2 {x,y,z}, t]`? Looks like not. If the numbers are rational, the results are generally expected but interesting. Either cycles of greater or lesser length are obtained, either $t/2$ occurs. But for reals I see some problems. For any initial data with machine precision we get "Half Total" sooner or later (here $t=1$): {{0.84474, 0.117199, 0.0380613}, {0.68948, 0.234397, 0.0761226}, {0.37896, 0.468795, 0.152245}, {0.24208, 0.0624103, 0.69551}, {0.48416, 0.124821, 0.39102}, {0.0316802, 0.750359, 0.217961}, {0.0633604, 0.500718, 0.435922}, {0.126721, 0.00143543, 0.871844}, {0.253441, 0.00287087, 0.743688}, {0.506883, 0.00574174, 0.487375}, {0.0137658, 0.0114835, 0.974751}, {0.0275315, 0.0229669, 0.949502}, {0.0550631, 0.0459339, 0.899003}, {0.110126, 0.0918678, 0.798006}, {0.220252, 0.183736, 0.596012}, {0.440505, 0.367471, 0.192024}, {0.118991, 0.265058, 0.615952}, {0.237981, 0.530115, 0.231903}, {0.475963, 0.060231, 0.463806}, {0.0480747, 0.879538, 0.0723872}, {0.0961494, 0.759076, 0.144774}, {0.192299, 0.518152, 0.289549}, {0.384598, 0.0363047, 0.579098}, {0.769195, 0.0726095, 0.158195}, {0.538391, 0.145219, 0.31639}, {0.0767815, 0.290438, 0.632781}, {0.153563, 0.580876, 0.265561}, {0.307126, 0.161751, 0.531122}, {0.614252, 0.323503, 0.062245}, {0.228504, 0.647006, 0.12449}, {0.457009, 0.294012, 0.24898}, {0.0859828, 0.411977, 0.50204}, {0.171966, 0.823954, 0.00408077}, {0.343931, 0.647907, 0.00816154}, {0.687862, 0.295815, 0.0163231}, {0.375725, 0.591629, 0.0326462}, {0.75145, 0.183258, 0.0652924}, {0.502899, 0.366516, 0.130585}, {0.00579834, 0.733032, 0.261169}, {0.0115967, 0.466064, 0.522339}, {0.0231934, 0.932129, 0.0446777}, {0.0463867, 0.864258, 0.0893555}, {0.0927734, 0.728516, 0.178711}, {0.185547, 0.457031, 0.357422}, {0.628906, 0.0859375, 0.285156}, {0.257813, 0.171875, 0.570313}, {0.515625, 0.34375, 0.140625}, {0.03125, 0.6875, 0.28125}, {0.0625, 0.375, 0.5625}, {0.125, 0.75, 0.125}, {0.25, 0.5, 0.25}, {{0.25, 0.5, 0.25}, "Half Total"}} and for higher precision we get "Total error" (sooner or later too): ClearAll[genNumsReals]; genNumsReals[tot_] := Block[{a = RandomReal[{0., tot}, WorkingPrecision -> 30], b}, b = RandomReal[{0., tot - a}, WorkingPrecision -> 30]; {a, b, tot - a - b}]; Clear[nextNums]; nextNums[nums_, tot_] := Which[ Unequal[Total@nums, tot], {nums, "Error Total"}, MemberQ[nums, tot/2], {nums, "Half Total"}, AnyTrue[nums, GreaterThan[tot/2]], SetPrecision[Mod[2 nums, tot], 30], True, SetPrecision[tot - 2 nums, 30]]; {{0.659587099030496875713494880795, 0.185853135754329321586581016710, 0.154559765215173802699924102495}, \ {0.319174198060993674630481109489, 0.371706271508658647739764546714, 0.309119530430347622118603112540}, \ {0.361651603878012650739037781022, 0.256587456982682704520470906573, 0.381760939139304755762793774920}, \ {0.276696792243974698521924437955, 0.486825086034634590959058186854, 0.236478121721390488474412450159}, \ {0.446606415512050602956151124090, 0.0263498279307308180818836262915, 0.527043756557219023051175099681}, \ {0.893212831024101205912302248180, 0.0526996558614616361637672525831, 0.0540875131144380461023501993623}, \ {0.786425662048202411824604496360, 0.105399311722923272327534505166, 0.108175026228876092204700398725}, \ {0.572851324096404823649208992720, 0.210798623445846544655069010332, 0.216350052457752184409400797449}, \ {0.145702648192809647298417985439, 0.421597246891693089310138020664, 0.432700104915504368818801594898}, \ {0.708594703614380705403164029121, 0.156805506216613821379723958671, 0.134599790168991262362396810204}, \ {{0.708594703614380705403164029121, 0.156805506216613821379723958671, 0.134599790168991262362396810204}, "Error Total"}} I think the problems are in rounding, but how to understand it?

For some reasons I need to consider the following iterative process.
We have three positive numbers with constant total such that $x+y+z == t$
Consider the iterative algorithm:

if the sum of numbers is not equal to $t$ - stop (with message "Total error")
if one of the numbers is $t/2$ - stop (with message "Half")
if one of the numbers is greater than $t/2$ (let it be $z$) - returns new numbers $2z-t,2x,2y$

if all numbers are less than $t/2$ - return the new numbers $t-2z,t-2x,t-2y$

ClearAll[genNumsReals]; genNumsReals[ tot_] :=
 Block[{a = RandomReal[{0., tot}], b},
  b = RandomReal[{0., tot - a}];
  {a, b, tot - a - b}];

ClearAll[genNumsRationals]; genNumsRationals[e_, tot_] :=
 Block[{a = Rationalize[RandomReal[{e, tot - e}], e], b},
  b = Rationalize[RandomReal[{e, tot - a - e}], e];
  {a, b, tot - a - b}];

Clear[nextNums]; nextNums[nums_, tot_] :=
 Which[
  Unequal[Total@nums, tot], {nums, "Error Total"},
  MemberQ[nums, tot/2], {nums, "Half Total"},
  AnyTrue[nums, GreaterThan[tot/2]], Mod[2 nums, tot],
  True, tot - 2 nums];

numsNest[init_, tot_, iterMax_ : 100] :=
  NestWhileList[nextNums[#, tot] &, init,
   (UnsameQ[##] && Length@Last[{##}] == 3) &, All, iterMax];

Firstly, is it equal to Mod[-2 {x,y,z}, t]? Looks like not.
If the numbers are rational, the results are generally expected but interesting.
Either cycles of greater or lesser length are obtained,
either $t/2$ occurs.

But for reals I see some problems. For any initial data with machine precision we get "Half Total" sooner or later (here $t=1$):

{{0.84474, 0.117199, 0.0380613}, {0.68948, 0.234397, 
  0.0761226}, {0.37896, 0.468795, 0.152245}, {0.24208, 0.0624103, 
  0.69551}, {0.48416, 0.124821, 0.39102}, {0.0316802, 0.750359, 
  0.217961}, {0.0633604, 0.500718, 0.435922}, {0.126721, 0.00143543, 
  0.871844}, {0.253441, 0.00287087, 0.743688}, {0.506883, 0.00574174, 
  0.487375}, {0.0137658, 0.0114835, 0.974751}, {0.0275315, 0.0229669, 
  0.949502}, {0.0550631, 0.0459339, 0.899003}, {0.110126, 0.0918678, 
  0.798006}, {0.220252, 0.183736, 0.596012}, {0.440505, 0.367471, 
  0.192024}, {0.118991, 0.265058, 0.615952}, {0.237981, 0.530115, 
  0.231903}, {0.475963, 0.060231, 0.463806}, {0.0480747, 0.879538, 
  0.0723872}, {0.0961494, 0.759076, 0.144774}, {0.192299, 0.518152, 
  0.289549}, {0.384598, 0.0363047, 0.579098}, {0.769195, 0.0726095, 
  0.158195}, {0.538391, 0.145219, 0.31639}, {0.0767815, 0.290438, 
  0.632781}, {0.153563, 0.580876, 0.265561}, {0.307126, 0.161751, 
  0.531122}, {0.614252, 0.323503, 0.062245}, {0.228504, 0.647006, 
  0.12449}, {0.457009, 0.294012, 0.24898}, {0.0859828, 0.411977, 
  0.50204}, {0.171966, 0.823954, 0.00408077}, {0.343931, 0.647907, 
  0.00816154}, {0.687862, 0.295815, 0.0163231}, {0.375725, 0.591629, 
  0.0326462}, {0.75145, 0.183258, 0.0652924}, {0.502899, 0.366516, 
  0.130585}, {0.00579834, 0.733032, 0.261169}, {0.0115967, 0.466064, 
  0.522339}, {0.0231934, 0.932129, 0.0446777}, {0.0463867, 0.864258, 
  0.0893555}, {0.0927734, 0.728516, 0.178711}, {0.185547, 0.457031, 
  0.357422}, {0.628906, 0.0859375, 0.285156}, {0.257813, 0.171875, 
  0.570313}, {0.515625, 0.34375, 0.140625}, {0.03125, 0.6875, 
  0.28125}, {0.0625, 0.375, 0.5625}, {0.125, 0.75, 0.125}, {0.25, 0.5,
   0.25}, {{0.25, 0.5, 0.25}, "Half Total"}}

and for higher precision we get "Total error" (sooner or later too):

ClearAll[genNumsReals]; genNumsReals[tot_] :=
 Block[{a = RandomReal[{0., tot}, WorkingPrecision -> 30], b},
  b = RandomReal[{0., tot - a}, WorkingPrecision -> 30];
  {a, b, tot - a - b}];

Clear[nextNums]; nextNums[nums_, tot_] :=
 Which[
  Unequal[Total@nums, tot], {nums, "Error Total"},
  MemberQ[nums, tot/2], {nums, "Half Total"},
  AnyTrue[nums, GreaterThan[tot/2]],
  SetPrecision[Mod[2 nums, tot], 30],
  True, SetPrecision[tot - 2 nums, 30]];

{{0.659587099030496875713494880795, 0.185853135754329321586581016710,
  0.154559765215173802699924102495}, \
{0.319174198060993674630481109489, 0.371706271508658647739764546714,
  0.309119530430347622118603112540}, \
{0.361651603878012650739037781022, 0.256587456982682704520470906573,
  0.381760939139304755762793774920}, \
{0.276696792243974698521924437955, 0.486825086034634590959058186854,
  0.236478121721390488474412450159}, \
{0.446606415512050602956151124090, 0.0263498279307308180818836262915,
  0.527043756557219023051175099681}, \
{0.893212831024101205912302248180, 0.0526996558614616361637672525831,
  0.0540875131144380461023501993623}, \
{0.786425662048202411824604496360, 0.105399311722923272327534505166,
  0.108175026228876092204700398725}, \
{0.572851324096404823649208992720, 0.210798623445846544655069010332,
  0.216350052457752184409400797449}, \
{0.145702648192809647298417985439, 0.421597246891693089310138020664,
  0.432700104915504368818801594898}, \
{0.708594703614380705403164029121, 0.156805506216613821379723958671,
  0.134599790168991262362396810204}, \
{{0.708594703614380705403164029121, 0.156805506216613821379723958671,
   0.134599790168991262362396810204}, "Error Total"}}

I think the problems are in rounding, but how to understand it?

POSTED BY: Denis Ivanov

3 Replies

Sort By:

Michael Rogers

Michael Rogers, Emory University

Posted 18 hours ago

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. 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'=1+2\delta$ 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.1102210^-16, -2.2204510^-16, 4.4408910^-16, 8.8817810^-16, -1.7763610^-15, -3.5527110^-15, -7.1054310^-15, -1.4210910^-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.5511210^-17 ) The error after the first step is twice this $-2^{53} = 1.11\times10^{-16}$. Example: Case $\delta = 0$. Remove the roundoff from an initial value less than 0.25 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

POSTED BY: Michael Rogers

Denis Ivanov

Posted 17 hours 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

Denis Ivanov

Posted 4 hours ago

So `SetPrecision[tot, Infinity]` helps a lot

POSTED BY: Denis Ivanov

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback

Problems with iterative process

Floating-point analysis

Example: Case $\delta \ne 0$

Example: Case $\delta = 0$.