# How to find an estimated root efficiently with binary stochastic estimation of the function?

I am trying to construct certain nice geometric patterns, and one particular feature I'm trying to optimize is the balance of black and white corresponding to point being inside, or outside a geometric region. The problem is that area computations over these parameterized regions can easily get too complicated to handle in any meaningful amount of time.

At the same time, membership queries are still rather straight-forward, and functions of area are continuous, so finding local roots per se shouldn't be awfully hard.

After the release of v10.2, I've got inspired by RandomPoint function to work around problems of computing these areas. So, I make plots like in this greatly simplified example:

Needs["ErrorBarPlots"];

Module[{samplesize, sample},
samplesize = 1000;
sample[x_] := RegionMember[Disk[{x, 0}], RandomPoint@Disk[]];
ErrorListPlot[
ParallelTable[{{x, # - 1/2},
ErrorBar[2 Sqrt[(# (1 - #))/samplesize]]} &[
Mean@Boole@Table[sample[x], {samplesize}]], {x, 0, 2, 1/20}],
Joined -> True]]


I plot the sample means (minus the target value 1/2, bringing roots to the axis) and two-sigma standard deviations of sample means, and visually inspect the graph. Then I pick a subsection, and increase sample size to sufficiently reduce error bars.

The question is: how to mechanize this process, and preferably take intelligently adaptive amount of samples, instead of just pulling these numbers out of a Stetson-Harrison hat? Essentially I want a FindRoot counterpart which could work with a function like sample (which returns either True or False) above, target mean, an initial starting point or interval, and some sort of an accuracy goal.

• This problem is sometimes called "stochastic root finding" or "stochastic approximation". Things like Newton's method will fail spectacularly on problems like these once the noise is high enough. I'm aware of two methods to deal with this: probabilistic bisection (nice introduction in this thesis) and the Robbins-Monro algorithm. This latter has an awful Wikipedia page, avoid it, except for references. This is more useful. Commented Jul 19, 2015 at 13:13
• Regarding Robbins-Monro, one thing that's not clear from the Wikipedia page is that there are many ways to do this, and how you do it precisely will have a very significant effect on convergence speed. Drop me a mail for a sloppy Mathematica implementation of probabilistic bisection ... Commented Jul 19, 2015 at 13:14
• @Szabolcs Thanks for these references! I'm pretty certain I've ran into them earlier, but happily forgot of them. In this case noise is of course very specific kind - not Gaussian noise, for instance. Commented Jul 19, 2015 at 13:18
• These might be off base but here are a couple of possibilities. (1) Hit (f[x]-k)^2 with findMinimum, letting k vary between -small and small. (Or use FindRoot directly on f[x]-k). Likely to do poorly in absence of decent starting point. (2) Select a random point p and direction, evaluate at p and along a ray in the selected direction. Once you get both positive and negative evalautions, use bisection to home in. (Maybe this is what you've been doing though.) Commented Jul 19, 2015 at 13:40
• @DanielLichtblau I'm a bit lost in the regard of f in your suggestions. Do you mean the a binary-valued (0 or 1) stochastic function similar to sample in my question? Commented Jul 19, 2015 at 13:45

This is a hack to create a function which interoperates at least on some level with FindRoot. Still, it's very much only a starting point. Primary goal of regularize is to attempt keeping values below and above sought root not jumping on the other side.

ClearAll@regularize;

regularize[f_Function, rootval_, minsigmas_, minsamples_Integer,
maxsamples_Integer, minsamesamples_Integer] :=
Module[{sigmaTest, regularized},
sigmaTest[succ_, total_] :=
(* maximum number of samples *)
total < maxsamples &&
(* minimum number of samples *)
(total < minsamples ||
(* number of samples required if all samples are the same *)
(total < minsamesamples && (succ == 0 || succ == total)) ||
(* minimum distance from root value in sigmas *)
Abs[succ/total - rootval] <
minsigmas Sqrt[(succ/total) ((total - succ)/total)/total]);
regularized[x_?NumericQ] := #1/#2 & @@
NestWhile[Apply[{#1 + Boole@f[x], #2 + 1} &], {0., 0.},
Apply@sigmaTest];
regularized]

With[{regularized =
regularize[RegionMember[Disk[{#, 0}], RandomPoint@Disk[]] &, 0.5,
6, 100, 500000, 100]},
FindRoot[regularized[x] == 0.5, {x, 0.5, 1.}, AccuracyGoal -> 3,
PrecisionGoal -> 3, EvaluationMonitor :> Print@x]]


0.5

1.

0.848923

0.808064

0.807856

0.807743

0.653872

0.806863

(* {x -> 0.806863} *)

True value of this root is about 0.80794551.

I don't like the amount of parameters and their somewhat arbitrary nature. I'm more concerned, though, of FindRoot treating these values as exact values, instead of probabilistic samples. Ideal solution on my mind would take accumulated probabilities into account when targeting its hunt of the highest-likelihood zero crossing, but I don't have energy to study that at this point.

This is my attempt at Probabilistic Bisection Algorithm, as described in the thesis linked by @Szabolcs in the comments.

ClearAll[pbaFindRoot];

Module[{initialDistribution, modifyDistribution, outputDistribution,
bernoulliTestPowerOneFilter},
initialDistribution[rangemin_, rangemax_] :=
N@{(rangemin + rangemax)/
2, {{{0, rangemin}, {(rangemin + rangemax)/2, 1/2}}, {{rangemax,
1/2}, {rangemax, 0}}}};

modifyDistribution[{_, lists_List}, p_] :=
Module[{list, sum0, sum1, len, median},
list = {
Sequence @@ ({1, 2 p} # & /@ #1),
Sequence @@ ({1, 2 (1 - p)} # & /@ #2)} & @@ lists;
sum1 = 0;
len = LengthWhile[list, (sum1 += Last@#; sum1 < 1/2) &] + 1;
sum0 = sum1 - list[[len, 2]];
median = Rescale[1/2, {sum0, sum1}, list[[len - 1 ;; len, 1]]];
{median,
{Append[list[[;; len - 1]], {median, 1/2 - sum0}],
Join[{list[[len]] - {0, 1/2 - sum0}},
list[[len + 1 ;;]]]}}];

outputDistribution[dist_List] :=
"Histogram", {#2[[2 ;; -2]]/Differences@Most@#1, Most@#1} & @@
Transpose[Flatten[Last@dist, 1]], 1, Infinity];

bernoulliTestPowerOneFilter[boolsample_Function, p_, gamma_: 0.5,
maxiterations_Integer: 1000] :=
#1/#2 > p & @@
NestWhile[Apply[{#1 + Boole@boolsample[], #2 + 1} &], {0, 0},
Apply[Function[{succ, n},
Abs[succ - n p] <
Sqrt[n (Log[n + 1] - Log[2] - Log[gamma])/
2]]], {2, 1}, maxiterations];

pbaFindRoot[boolsample_Function, targetp_, {rangemin_, rangemax_},
iter_: 100, distskew_: 0.1, gamma_: 0.5, poweronemaxiter_: 1000] :=
outputDistribution@Nest[
With[{v = First@#},
modifyDistribution[#,
0.5 -
distskew + (2 distskew) Boole[
bernoulliTestPowerOneFilter[boolsample[v] &, targetp,
gamma, poweronemaxiter]]]] &,
initialDistribution[rangemin, rangemax], iter]];

ClearAll[dist];

dist = pbaFindRoot[
Not@RegionMember[Disk[{#, 0}], RandomPoint@Disk[]] &, 0.5, {0, 2},
200, 0.1, 0.5, 5000];

{DiscretePlot[PDF[dist, x], {x, 0, 2, .001}, PlotRange -> All],
Median@dist, StandardDeviation@dist}


It is noteworthy that this implementation doesn't return a plain value, but actual distribution generated by the PBA algorithm.

Well, it works. Sort of. It would seem at least with some monotonically increasing functions with the root (targetp) in the range. The most opaque part of this code is probably modifyDistribution` which maintains lists of distribution data in a way which doesn't destroy performance right away.