Rounding Randomly - Reasonably Right?

2025-08-21

Introduction

If you have access to a terminal right now, do me a favor - open it up and type in the following command:

python3 -c "print(0.1 + 0.2)"

If you are surprised by the output, don't worry. In fact, your python interpreter is working exactly as it should, as specified by the IEEE 754 standard for floating-point arithmetic. Programmers sometime assume that floating-point numbers are infinitely precise, but that is simply not the case. To quote Adam Sawicki's blog on the topic (Myths About Floating-Point Numbers)[0]:

They are not exact

It is not true that 2.0 + 2.0 can give 3.99999. It will always be 4.0. They are exact to the extent of their limited range and precision.

The following article does not go deep into explaining why this is the case and how exactly floating-point arithmetic works, but if you are interested in that I definitely recommend the post mentioned above, as well as this very useful primer: What Is Floating-Point Arithmetic[1].

Instead of that, this article is about an alternative way of handling the rounding mechanism - through a method called stochastic rounding. I was first introduced to this method during my undergrad studies and was fascinated by it. Two years later, I thought it would be a good idea to share that fascination with the rest of the world along with some of my views on it. And so in this article we will be looking into the implementation of this method in c, discussing if it's even any good or not, and if so, what applications it might prove useful in. I hope you enjoy :).

Stochastic Rounding

Quoting from Croci et al. (2021)[3]:

Stochastic rounding randomly maps a real number x to one of the nearest values in a finite precision number system. The probability of choosing either of these two numbers is 1 minus their relative distance to x.

Put more simply, it is a non-deterministic rounding method which can help to cancel out rounding errors over many rounding operations. In a more formal and mathematical sense, stochastic rounding can be defined as follows:

where \(x\) is some real value, \(\text{RZ}(x)\) is the rounding mode round-toward-zero, \(\text{RA}(x)\) is round-away-from-zero rounding mode, and \(P \in [0, 1]\) is a random number from a uniform distribution.

Note that this particular flavor of stochastic rounding is a little different from another implementation: the other form rounds the value up or down with equal probability(1/2), but this form takes into account the distance between the value and its up or down representations.

Implementation of binary32 stochastic rounding

NB: The c file I used for writing and testing the implementation is available here.

The following function takes as input a double precision floating point number and returns a single precision floating point number after applying stochastic rounding:

float nearestup(double val) {
	float nearest = (float)val;
	if (nearest > val){
		return nearest;
	}
	return nextafterf(nearest, INFINITY);
}

float nearestdown(double val) {
	float nearest = (float)val;
	if (nearest < val){
		return nearest;
	}
	return nextafterf(nearest, -INFINITY);
}

float srounding(double val) {
	if (val == (float)val) {
		return (float)val;
	}
	float down = nearestdown(val);
	float up = nearestup(val);
	double randprob = (double)rand()/(double)RAND_MAX;

	if (randprob < ((val - down)/(up - down))) {
		return up;
	} else {
		return down;
	}
}

The important built-in C functions that were used in the implementation of this function were rand() which returns a random number between 0 and RAND_MAX and this was used to generate the random number between 0 and 1. The nextafterf() function was also used to get the smallest integer greater than or equal to x and the largest integer less than or equal to x.

Usefulness

Stochastic Rounding forces the rounding errors to be random variables with zero mean.

The following proof is taken from Connolly et al.(2020) [4]

Recall the definition of SR:

Since \(P\) is a random variable, \(SR(x)\) and \(\delta = x - \text{SR}(x)\) are also random variables, the expected value of \(SR(x)\) would be:

Which results in:

And so:

The expected value of the rounding error is zero.

The benefit of stochastic rounding is that it maintains some information that would otherwise be lost in deterministic rounding methods. Verified in many studies, SR is shown to have a better error bound than its more popular alternate Round to Nearest(RN) while computing the inner product of two vectors, multiplying matrices and vectors, as well as summation.

Testing π

We will now test our implementation by rounding pi (as represented in double precision) using stochastic rounding to binary32 precision.

int main(int argc, char** argv){
  double sample = M_PI;
  const long int tries = 5000000;

  double rounded_val = 0;
  float rounded = 0;
  for (int i=1; i<=tries; i++) {
    rounded = stocrounding(sample);
    rounded_val += rounded;
  }
  rounded_val = rounded_val/tries;

  printf("Value being rounded in binary64:                         %.60f \n", sample);
  printf("Represenation in binary32:                               %.60f \n", (float)sample);
  printf("Average of all rounded values using stochastic rounding: %.60f \n", rounded_val);
}

The output of the above code is as follows:

Value being rounded in binary64:                         3.141592653589793115997963468544185161590576171875000000000000
Represenation in binary32:                               3.141592741012573242187500000000000000000000000000000000000000
Average of all rounded values using stochastic rounding: 3.141592653529119427702198663610033690929412841796875000000000

Usefulness

Imagine a large city, with a network of advanced temperature sensors spread out all over. Each sensor is capable of measuring the ambient temperature up to a 64 bit floating-point value, but there is a bottleneck however: the network connecting these sensors to the central data hub can only handle smaller 32-bit numbers. This forces each sensor to simplify, or round, its precise measurement before sending it.

The common approach would be Round-to-Nearest. Let's see how this might affect measurements through an example: three sensors measure temperatures of 20.55 °C, 21.55 °C, and 22.55 °C. The true average of these is exactly 21.55 °C.

If the 32-bit format can only store one decimal place, a "round-half-up" rule would round all these values up. The sensors would send 20.6 °C, 21.6 °C, and 22.6 °C, which would be averaged to 21.6 °C. The bias accumulates, and there is a drift from the true values.

This is where SR might prove to be useful. For a value like 20.55 which is 55% of the way between 20.5 and 20.6, the algorithm gives it a 55% chance of being rounded up to 20.6 and a 45% chance of being rounded down. What this means is that if the value was rounded 100 times, about 55 of those values would 20.6 and about 45 would be 20.5. Of course, rounding 100 times might simply not be feasible, but perhaps rounding 3 or so times would be: In that case, the sensors could return the following array of values: - [20.6, 20.5, 20.6], Average: 20.566 - [21.6, 21.5, 21.6], Average: 21.566 - [22.6, 22.5, 22.6], Average: 22.566 With the overall average being: 21.566.

This is obviously quite a stupid example, but it is helpful in understanding the tasks where SR might prove to be very useful. Let's move to a real computational scenario next.

Gradient Descent

To probe whether stochastic rounding (SR) helps in a realistic-looking training loop, I ran a small gradient-descent microbenchmark. The goal is not to claim production-level speedups, but to show whether SR can reduce the accumulated rounding bias that arises from many tiny weight updates.

The code:

void gradient_descent_test() {
  const int epochs = 1000;
  const double learning_rate = 1e-4;
  const int num_weights = 1000;
  const int trials = 20;

  double *grads = malloc((size_t)epochs * num_weights * sizeof(double));
  double *initial_weights = malloc(num_weights * sizeof(double));
  double *weights_true = malloc(num_weights * sizeof(double));
  float *weights_std = malloc(num_weights * sizeof(float));

  for (int w = 0; w < num_weights; w++) {
    initial_weights[w] = 0.1 * (rand01() - 0.5); // [-0.05, 0.05]
    weights_true[w] = initial_weights[w];
    weights_std[w] = (float)weights_true[w];
  }

  for (int epoch = 0; epoch < epochs; epoch++) {
    for (int w = 0; w < num_weights; w++) {
      double base_grad = 1e-6 * sin(epoch * 0.01 + w * 0.1);
      double noise = 5e-7 * (rand01() - 0.5);
      grads[epoch * num_weights + w] = base_grad + noise;
    }
  }

  for (int epoch = 0; epoch < epochs; epoch++) {
    for (int w = 0; w < num_weights; w++) {
      double gradient = grads[epoch * num_weights + w];
      weights_true[w] += learning_rate * gradient;
      double update = learning_rate * gradient;
      weights_std[w] += (float)update;
    }
  }

  double total_error_std = 0.0;
  double total_error_sr = 0.0;

  for (int trial = 0; trial < trials; trial++) {
    float *weights_sr = malloc(num_weights * sizeof(float));

    for (int w = 0; w < num_weights; w++) {
      weights_sr[w] = (float)initial_weights[w]; // Same as standard
    }

    for (int epoch = 0; epoch < epochs; epoch++) {
      for (int w = 0; w < num_weights; w++) {
        double gradient = grads[epoch * num_weights + w];
        double update = learning_rate * gradient;
        weights_sr[w] = srounding((double)weights_sr[w] + update);
      }
    }

    for (int w = 0; w < num_weights; w++) {
      double error_sr = fabs(weights_sr[w] - weights_true[w]);
      total_error_sr += error_sr * error_sr;
    }

    free(weights_sr);
  }

  for (int w = 0; w < num_weights; w++) {
    double error_std = fabs(weights_std[w] - weights_true[w]);
    total_error_std += error_std * error_std;
  }

  double rmse_std = sqrt(total_error_std / num_weights);
  double rmse_sr = sqrt(total_error_sr / (num_weights * trials));

  printf("Training simulation: %d epochs, %d weights\n", epochs, num_weights);
  printf("Learning rate: %.1e\n", learning_rate);
  printf("Gradient scale: ~1e-6 to 1e-7\n");
  printf("\n");
  printf("Standard rounding RMSE:   %.2e\n", rmse_std);
  printf("Stochastic rounding RMSE: %.2e\n", rmse_sr);

  if (rmse_sr < rmse_std) {
    printf("Improvement factor:       %.2fx\n", rmse_std / rmse_sr);
  } else {
    printf("Standard is better by:    %.2fx\n", rmse_sr / rmse_std);
  }

  free(weights_true);
  free(weights_std);
  free(initial_weights);
  free(grads);
  printf("\n");
}

To keep the comparison fair, I pre-generated the entire sequence of per-step gradients and reused that exact gradient stream for the double-precision reference and for every SR trial; this way any difference between runs comes only from rounding behavior, not from different gradients or initial conditions. The reference (weights_true) was updated in double precision; the deterministic baseline (weights_std) kept a single-precision copy updated with round-to-nearest casts; the SR runs kept the weights in single precision but stochastically rounded the result of each accumulation step using srounding((double)weights_sr + update). I repeated the SR experiment across independent trials and computed per-weight RMSE versus the double reference.

These were the results I achieved:

=== Gradient Descent Test ===
Training simulation: 1000 epochs, 1000 weights
Learning rate: 1.0e-04
Gradient scale: ~1e-6 to 1e-7

Standard rounding RMSE:   1.32e-08
Stochastic rounding RMSE: 1.16e-08
Improvement factor:       1.14x

The absolute errors are very small, so the relative improvement is modest but measurable. Why does SR help here? SR makes each local rounding error unbiased in expectation; when an algorithm performs many such roundings, those zero-mean local errors tend to cancel on average, whereas deterministic rounding can introduce a systematic bias that accumulates. The result is a smaller expected deviation from the high-precision trajectory. Because SR is stochastic, however, single trials can be better or worse than round-to-nearest: the improvement emerges in the average across trials, so one must report both the mean and the variability.

References

  1. Sawicki, Adam (2021). Myths about Floating-Point Numbers. Blog post

  2. Higham, N. J. (2020). What is Floating-Point Arithmetic? Blog post

  3. Higham, N. J. (2020). What is Stochastic Rounding? Blog post

  4. Croci, M., Fasi, M., Higham, N. J., Mary, T., & Mikaitis, M. (2022). Stochastic rounding: Implementation, error analysis and applications. Royal Society Open Science, 9(211631). DOI

  5. Connolly, M. P., Higham, N. J., & Mary, T. (2020). Stochastic rounding and its probabilistic backward error analysis (MIMS EPrint: 2020.12).