<p>Let <span class="math-container">$x_1, \ldots, x_n \in \mathbb{R}$</span>, and let <span class="math-container">$y = f(x_i)$</span> be an arithmetic circuit in the <span class="math-container">$x_i$</span>'s.  That is, <span class="math-container">$f$</span> is a circuit of negate, add, subtract, and multiply gates.  Let <span class="math-container">$X_i$</span> be floating point random variables with <span class="math-container">$E(X_i) = x_i$</span>, and let <span class="math-container">$Y$</span> be the result of evaluating <span class="math-container">$f(X_i)$</span> using <a href="https://en.wikipedia.org/wiki/Rounding#Stochastic_rounding" rel="nofollow noreferrer">stochastic rounding</a> to the two nearest floating point numbers at every step.  That is, at each step we evaluate <span class="math-container">$\operatorname{op}(a,b)$</span> in infinite precision, then randomly choose one of the two neighboring floating point values so that the mean choice is the exact answer.</p>
<p><strong>Informal question:</strong> If we sample many times from <span class="math-container">$Y$</span>, can we get a PAC estimate of <span class="math-container">$y$</span> from the sample mean and variance?</p>
<p>Here's why I think this might be doable, at least in some cases:</p>
<ol>
<li><p><a href="https://epubs.siam.org/doi/pdf/10.1137/20M1334796?download=true" rel="nofollow noreferrer">Connolly et al. 2021, Stochastic rounding and its probabilistic backward error analysis</a> show that if the circuit <span class="math-container">$f$</span> is a tree, the errors at each step are mean independent (mean zero conditional on previous errors).  They use this to derive backward error estimates that are roughly as good as assuming the errors are entirely independent.</p>
</li>
<li><p>If <span class="math-container">$f$</span> is linear, then <span class="math-container">$E(Y) = y$</span>.</p>
</li>
<li><p>The first problem where mean independence breaks is squaring: <span class="math-container">$E(X^2) \ne E(X)^2$</span>.  However, if <span class="math-container">$X = x(1 + E)$</span> with <span class="math-container">$|E| \le \epsilon$</span> for <span class="math-container">$\epsilon \ll 0$</span>, then
<span class="math-container">$$X^2 = x^2(1+2E+E^2) = x^2(1+2E) + O(\epsilon^2)$$</span>
Ideally, this would mean that if the circuit <span class="math-container">$f$</span> has low multiplication depth, the linear error term dominates the distribution of <span class="math-container">$Y$</span>, and that term is mean zero.</p>
</li>
</ol>
<p>The circuit I'm interested in is not particularly stable, so I don't expect to be able to prove anything about its error properties.  So the ideal would be to be able to read off an accuracy estimate from samples from <span class="math-container">$Y$</span>.  Unfortunately, the &quot;non-mean-zero term is small&quot; intuition from (3) might be impossible to detect from the distribution of <span class="math-container">$Y$</span>.</p>
<p><strong>More precise question:</strong> Say we sample <span class="math-container">$y_1, \ldots, y_n \sim Y$</span>, and fix <span class="math-container">$\lambda &gt; 1$</span>.  Can we bound the probability <span class="math-container">$p$</span> s.t. <span class="math-container">$y \in [m-\lambda s, m+\lambda s]$</span>, where <span class="math-container">$m = \frac{1}{n}\sum_i y_i, s = \frac{1}{n}\sum_i (y_i-m)^2$</span> are the sample mean and variance?</p>
<p><strong>Floating point details:</strong> I'm happy to ignore exponent overflow/underflow, so we can take as our representable floating point values <span class="math-container">$\mathbb{F} = \{0\} \cup \pm\{2^{p-1}, \ldots, 2^p-1\} 2^\mathbb{Z}$</span> for some bit precision <span class="math-container">$p \in \mathbb{N}$</span> (<span class="math-container">$p = 53$</span> for double precision).  We can assume oracle access to <span class="math-container">$X_i$</span>, either because <span class="math-container">$x_i \in \mathbb{F}$</span> or because we are given a routine to sample each <span class="math-container">$X_i$</span>.</p>

CC BY-SA 4.0
2022-02-03T14:58:16.753Z
Maybe [this chaotic quadratic](https://en.wikipedia.org/wiki/Logistic_map) would be a good candidate for a counter-example?  E.g. have one input $x_1$, define $f_0(x_1)=x_1$ and define $f_n(z) = 4 f_{n-1}(z)(1-f_{n-1}(z))$.  This is chaotic so for large enough $n$ you will have $|f_n(1/3+\epsilon) - f_n(1/3-\epsilon)|$ being large even for very small $\epsilon$.  Evaluating in python on my machine with $n=200$ I see $f_n(2^{-50}) \approx 0.16722033851971366$, while $f_n(2^{-51}) \approx 0.8150761827512897$.  And of course $f_n(0)=0$.
Neal Young
CC BY-SA 4.0
2022-02-03T23:47:52.840Z
@NealYoung: But wouldn't that produce large variance in $Y$?  The question is roughly whether, if $\operatorname{var}(Y)$ is small, one can conclude we have a good estimate for $y$.
Geoffrey Irving
CC BY-SA 4.0
2022-02-04T04:39:14.863Z
Oh, yes, I see what you mean.  A bad example would have to be tightly concentrated, but around the wrong answer.  Try iterating the map $f(x) = -0.23 x^3 - 0.54 x^2 + 1.77x$.  It should have $f(0) = 0$ and $f(1) = 1$, and iterating starting at $x=0.001$ or $x=-0.001$, say (very near 0) should converge to 1 instead of 0.  So probably this can be converted into a counterexample (by scaling so that the fixed point 0 is moved to, say, 1/3, which is not representable, then each sample of $f_n(1/3)$ (the $n$th iteration of $f$, for large but fixed $n$) should be near 1, even though $f_n(1/3) = 1/3$.
Neal Young
CC BY-SA 4.0
2022-02-04T04:42:11.737Z
For the lazy, [here](https://www.google.com/search?q=plot+-0.23x%5E3-0.54x%5E2%2B1.77x%2C+x&client=safari&rls=en&sxsrf=APq-WBsGQvPVabPizM7KC7Uxtix8mvfWyA%3A1643948735303&ei=v6r8YfaEEqaCwbkP0oyPgAc&ved=0ahUKEwi24JWjmuX1AhUmQTABHVLGA3AQ4dUDCA0&uact=5&oq=plot+-0.23x%5E3-0.54x%5E2%2B1.77x%2C+x&gs_lcp=Cgdnd3Mtd2l6EAM6BwgAEEcQsANKBAhBGABKBAhGGABQqAtY-lBgilhoAXACeACAAeIBiAGsB5IBBTYuMi4xmAEAoAEByAEIwAEB&sclient=gws-wiz) is a plot of $f(x) = -0.23x^3 - 0.54x^2+1.77x$.
Neal Young
CC BY-SA 4.0
2022-02-04T09:42:05.703Z
@NealYoung: That's a nice idea!  It does seem like that sinks the hope of reading off accuracy from just the final answer.  In that example the computation would still be backward stable, but plausibly one could a variant with a slightly higher degree polynomial that doesn't get backward stability either.
Geoffrey Irving
<randomized-algorithms>
Evaluating arithmetic circuits with stochastic rounding
Geoffrey Irving
# Evaluating arithmetic circuits with stochastic rounding

(Asked by: username_0 on 2022-02-01)

Let $x\_1, \ldots, x\_n \in \mathbb{R}$, and let $y = f(x\_i)$ be an arithmetic circuit in the $x\_i$'s. That is, $f$ is a circuit of negate, add, subtract, and multiply gates. Let $X\_i$ be floating point random variables with $E(X\_i) = x\_i$, and let $Y$ be the result of evaluating $f(X\_i)$ using [stochastic rounding](https://en.wikipedia.org/wiki/Rounding#Stochastic_rounding) to the two nearest floating point numbers at every step. That is, at each step we evaluate $\operatorname{op}(a,b)$ in infinite precision, then randomly choose one of the two neighboring floating point values so that the mean choice is the exact answer.

**Informal question:** If we sample many times from $Y$, can we get a PAC estimate of $y$ from the sample mean and variance?

Here's why I think this might be doable, at least in some cases:

1. [Connolly et al. 2021, Stochastic rounding and its probabilistic backward error analysis](https://epubs.siam.org/doi/pdf/10.1137/20M1334796?download=true) show that if the circuit $f$ is a tree, the errors at each step are mean independent (mean zero conditional on previous errors). They use this to derive backward error estimates that are roughly as good as assuming the errors are entirely independent.
2. If $f$ is linear, then $E(Y) = y$.
3. The first problem where mean independence breaks is squaring: $E(X^2) \ne E(X)^2$. However, if $X = x(1 + E)$ with $|E| \le \epsilon$ for $\epsilon \ll 0$, then
   $$X^2 = x^2(1+2E+E^2) = x^2(1+2E) + O(\epsilon^2)$$
   Ideally, this would mean that if the circuit $f$ has low multiplication depth, the linear error term dominates the distribution of $Y$, and that term is mean zero.

The circuit I'm interested in is not particularly stable, so I don't expect to be able to prove anything about its error properties. So the ideal would be to be able to read off an accuracy estimate from samples from $Y$. Unfortunately, the "non-mean-zero term is small" intuition from (3) might be impossible to detect from the distribution of $Y$.

**More precise question:** Say we sample $y\_1, \ldots, y\_n \sim Y$, and fix $\lambda > 1$. Can we bound the probability $p$ s.t. $y \in [m-\lambda s, m+\lambda s]$, where $m = \frac{1}{n}\sum\_i y\_i, s = \frac{1}{n}\sum\_i (y\_i-m)^2$ are the sample mean and variance?

**Floating point details:** I'm happy to ignore exponent overflow/underflow, so we can take as our representable floating point values $\mathbb{F} = \{0\} \cup \pm\{2^{p-1}, \ldots, 2^p-1\} 2^\mathbb{Z}$ for some bit precision $p \in \mathbb{N}$ ($p = 53$ for double precision). We can assume oracle access to $X\_i$, either because $x\_i \in \mathbb{F}$ or because we are given a routine to sample each $X\_i$.

**Tags:** randomized-algorithms

> **Comments:**
> **username_1**: Maybe [this chaotic quadratic](https://en.wikipedia.org/wiki/Logistic_map) would be a good candidate for a counter-example?  E.g. have one input $x_1$, define $f_0(x_1)=x_1$ and define $f_n(z) = 4 f_{n-1}(z)(1-f_{n-1}(z))$.  This is chaotic so for large enough $n$ you will have $|f_n(1/3+\epsilon) - f_n(1/3-\epsilon)|$ being large even for very small $\epsilon$.  Evaluating in python on my machine with $n=200$ I see $f_n(2^{-50}) \approx 0.16722033851971366$, while $f_n(2^{-51}) \approx 0.8150761827512897$.  And of course $f_n(0)=0$.  
> **username_0**: @NealYoung: But wouldn't that produce large variance in $Y$?  The question is roughly whether, if $\operatorname{var}(Y)$ is small, one can conclude we have a good estimate for $y$.  
> **username_1**: Oh, yes, I see what you mean.  A bad example would have to be tightly concentrated, but around the wrong answer.  Try iterating the map $f(x) = -0.23 x^3 - 0.54 x^2 + 1.77x$.  It should have $f(0) = 0$ and $f(1) = 1$, and iterating starting at $x=0.001$ or $x=-0.001$, say (very near 0) should converge to 1 instead of 0.  So probably this can be converted into a counterexample (by scaling so that the fixed point 0 is moved to, say, 1/3, which is not representable, then each sample of $f_n(1/3)$ (the $n$th iteration of $f$, for large but fixed $n$) should be near 1, even though $f_n(1/3) = 1/3$.  
> **username_1**: For the lazy, [here](https://www.google.com/search?q=plot+-0.23x%5E3-0.54x%5E2%2B1.77x%2C+x&client=safari&rls=en&sxsrf=APq-WBsGQvPVabPizM7KC7Uxtix8mvfWyA%3A1643948735303&ei=v6r8YfaEEqaCwbkP0oyPgAc&ved=0ahUKEwi24JWjmuX1AhUmQTABHVLGA3AQ4dUDCA0&uact=5&oq=plot+-0.23x%5E3-0.54x%5E2%2B1.77x%2C+x&gs_lcp=Cgdnd3Mtd2l6EAM6BwgAEEcQsANKBAhBGABKBAhGGABQqAtY-lBgilhoAXACeACAAeIBiAGsB5IBBTYuMi4xmAEAoAEByAEIwAEB&sclient=gws-wiz) is a plot of $f(x) = -0.23x^3 - 0.54x^2+1.77x$.  
> **username_0**: @NealYoung: That's a nice idea!  It does seem like that sinks the hope of reading off accuracy from just the final answer.  In that example the computation would still be backward stable, but plausibly one could a variant with a slightly higher degree polynomial that doesn't get backward stability either.  

---

## 0 answers

