Efficient Programmable Random Variate Generation Accelerator from Sensor Noise

by   James Timothy Meech, et al.
University of Cambridge

We introduce a method for non-uniform random number generation based on sampling a physical process in a controlled environment. We demonstrate one proof-of-concept implementation of the method that reduces the error of Monte Carlo integration of a univariate Gaussian by 1068 times while doubling the speed of the Monte Carlo simulation. We show that the supply voltage and temperature of the physical process must be controlled to prevent the mean and standard deviation of the random number generator from drifting.


page 1

page 2

page 3

page 4


A System for Generating Non-Uniform Random Variates using Graphene Field-Effect Transistors

We introduce a new method for hardware non-uniform random number generat...

On the Sampling Problem for Kernel Quadrature

The standard Kernel Quadrature method for numerical integration with ran...

Hamiltonian Monte Carlo for efficient Gaussian sampling: long and random steps

Hamiltonian Monte Carlo (HMC) is a Markov chain algorithm for sampling f...

A Scrambled Method of Moments

Quasi-Monte Carlo (qMC) methods are a powerful alternative to classical ...

Monte Carlo integration with adaptive variance reduction: an asymptotic analysis

The crude Monte Carlo approximates the integral S(f)=∫_a^b f(x) dx ...

Monte Carlo Convolution for Learning on Non-Uniformly Sampled Point Clouds

We propose an efficient and effective method to learn convolutions for n...

Fast random number generator based on optical physical unclonable functions

We propose an approach for fast random number generation based on homema...

I Introduction

Current software-based methods of non-uniform random variate generation are slow and inefficient [18][19][15][17]

. We present a programmable system capable of generating Gaussian random variates by extracting the noise properties of a MEMS sensor. We demonstrate its principle and application. Sampling a random physical process with a given distribution provides a continuous random variable with a theoretically unlimited sample rate. In a hardware implementation, the sample rate of the analog to digital convertor limits random number generation rate.

I-a Generating Uniform Random Variates Is Easy

Uniform random numbers are generated for cryptography [13]. Gaussian random variate generation is typically an order of magnitude slower and less efficient than uniform random variate generation [18]. Table I

compares several state-of-the-art methods of generating non-uniform random variates with a Gaussian distribution. We propose a method superior to all of the state-of-the-art methods in terms of sample rate and efficiency, consisting of a physical noise source and an analog to digital convertor.

Device Sample rate (Gsample/s) Efficiency (Msample/J)
CPU 0.89 3.17
GPU 12.90 107.52
MPPA 0.86 61.48
FPGA 12.10 403.20
ADC 14.00 645.10
TABLE I: Comparison of state-of-the-art methods for 16-bit Gaussian random variate generation [18]. The final row of the table represents seven dual-channel 1 Gsample/s ADS54J60s sampling each axis of five ADXL335s. Assuming that the samples move into the CPU memory in negligible time.

I-B Generating Non-Uniform Random Variates Is Harder

The inversion method and accept reject method are used in software for generating samples from non-uniform random variates [6]. Let and be uniform and non-uniform random variates respectively and

the analytical closed-form solution for the inverse cumulative distribution function 

[6]. Algorithm 1 shows the inversion method. The inversion method requires that the inverse cumulative distribution function has an analytical closed-form solution. The Gaussian distribution has no analytical closed-form solution for its inverse cumulative distribution function so cannot be used with the inversion method [17]. The accept reject method must be used instead. Let and be independent and the density on of . Where is -dimensional Euclidean space. Algorithm 2 shows the accept reject method. The accept reject method requires more math operations than the inversion method to transform a sample [6]. This causes it to take more clock cycles to compute and therefore, more time to transform each sample. When using the accept reject method, samples deviating from the desired distribution are rejected [6].

Result: Sample from non-uniform random variate
Generate uniform [0,1] random variate RETURN
Algorithm 1 Inversion method.
initialization initialization Result: Sample from non-uniform random variate
while  do
       Generate uniform [0,1] random variate Generate uniform [0,1] random variate Set RETURN
end while
Algorithm 2 Accept reject method.

I-C Uses Of Non-Uniform Random Variates

Non-uniform random variate generators are fundamental to applications employing Monte Carlo Methods [16][8][11], population balance modelling of the crystallisation process [1], generating phylogenetic trees [12], drug discovery [2], ray tracing [3], communication channel emulation [5], financial computing [20], local area networks [14], modelling of manufacture systems [10] and measuring healthcare strategy effectiveness [4].

When conducting Bayesian inference in probabilistic machine learning, we must evaluate Bayes’ theorem to calculate the probability that a belief

is true given new data , we denote this as . To do this we need the probability that the belief is true regardless of our data , the probability that the data is true given the belief and the probability that the data is true regardless of the belief . We will refer to as the prior, as the likelihood, as the marginal likelihood and as the posterior [9].


We calculate the marginal likelihood by integrating the joint density  [9]

. In practice, the analytical calculation of the marginal likelihood is impossible for all but the simplest joint distributions 

[9]. We instead sample from the joint distribution to obtain summary statistics that we can use to describe it [9]

. These random samples from bespoke probability distributions must be produced using a non-uniform random variate generator.

I-D Contributions

  1. The idea that physical noise sources such as MEMS sensors can be used as non-uniform random variate generators (Section I).

  2. Estimation of performance increase and error reduction achieved by using such a generator for Monte Carlo integration (Section II).

  3. Investigation of the impact of temperature and supply voltage on the noise distribution obtained from a commercial MEMS sensor (Section III-A).

Ii Motivating Example

We performed Monte Carlo integration of a Gaussian with a mean of 0 and standard deviation of 1 using samples from a Gaussian with the same parameters. We ran the experiment with a Gaussian generated by the C++ random library and repeated it with the random number generation rate adjusted to match that of our proposed non-uniform random variate generator based on physical noise sources. We performed the same integration using samples from a uniform distribution with a minimum of -3 and a maximum of 3. Let

be the error of the integration, be the time taken by the integration, be the number of random samples and be the distribution (either uniform or Gaussian). Let samps be the array of random numbers, be the area, be the rectangle base length, be the rectangle height and

the probability density function of the Gaussian for integration. Algorithm 

3 shows the integration scheme that we used. We repeated each integration process 1000 times and calculated the average. Figure 1 shows that the proposed hardware random number generator outperforms the C++ uniform random number generator after 1000 samples and reduces the error by after 1 million samples. The proposed hardware random number generator always performs the task at least twice as fast as the C++ Gaussian random number generator.

Result: Error and time
Timer start Generate random samples from distribution Sort random samples for All pairs of samples do
       b = samps[i] - samps[i-1] h = (f(samps[i]) + f(samps[i-1]))/2 A += b*h
end for
e = abs(1-A) Timer stop t = stop - start RETURN e,t
Algorithm 3 MC Integration.
Fig. 1: The error of Monte Carlo integration plateaus for uniform random numbers but decreases with a power relation for Gaussian random numbers.

Iii Methodology

A non-uniform random variate generator based on physical noise sources must have negligible drift of the mean and standard deviation over time. Drift would cause errors in calculations using the output of the random number generator. Any environmental parameter that causes non-negligible drift must be physically controlled, for example using a voltage regulator to control voltage and a Peltier device to control temperature. We investigated the dependence of 100,000 sample distributions on voltage and temperature. We sampled the z-axis of a MEMS accelerometer (we used the accelerometer in the Bosch BMX055) to obtain the distributions.

Iii-a Temperature-Controlled Experiments

Figure 2 shows the experimental setup. We placed the microcontroller, accelerometer, tilt and rotate stage and vibration isolation platform inside a Binder MK56 thermal chamber. We connected a microcontroller to the sensor via I2C for a 1154 Hz sample rate. We used a Keithly 2450 source measurement unit to power the sensor and measure the current drawn. We set the chamber temperature to 25 C and allowed 30 minutes for the temperature of the sensor to equilibrate whilst constantly sampling z-axis acceleration values from it. We then sampled 100,000 values from the BMX055 sensor at a 3.6 V supply voltage. We repeated this for all the voltages in the range of 3.6 to 1.4 V with a 0.2 V decrement. We then repeated this process for the temperature from 25 down to -5 C with a decrement of 5 C.

Fig. 2: Experimental setup, powering the sensor and the I2C interface separately allowed more accurate measurement of the power consumed by the sensor.

Iii-B Quantisation Investigation

We investigated the effect of quantisation on the Kullback–Leibler (KL) divergence between a discrete distribution and its ideal fitted curve. We used the MATLAB normrnd function to generate 100,000 values from a Gaussian distribution with the same mean and standard deviation as the BMX055 z-axis at 2.6 V and 10 C. We then discretized the values into various numbers of bins, fitted a Gaussian distribution to them and calculated the KL divergence between the fitted distribution and the actual distribution.

Iv Results and Discussion

We calculated the KL divergence between two discrete distributions using the following equation where and are discrete probability distributions, is a given sample value and the sample space [7]:


Figure 3 shows the BMX055 z-axis acceleration distribution at 2.6 V and 10 C. We found that the KL divergence between the distribution from the BMX055 z-axis and its fitted Gaussian (0.00263) was more than an order of magnitude smaller than the equivalent result for a MATLAB generated distribution (0.0392). We rounded the MATLAB generated floats for comparison with the BMX055 generated integers. The KL divergence between a MATLAB generated uniform distribution and its fitted Gaussian is 0.116 for reference. We averaged the KL divergence calculations over 100 distributions to account for the random variations in the measurement. The numbers generated by the sensor are closer to an ideal Gaussian distribution than those generated by MATLAB.

Fig. 3: Distribution of BMX055 accelerometer z-axis noise and closest fitted distribution. The KL divergence or difference between the data and the fitted distribution is 0.00263 demonstrating that the accelerometer noise closely matches a Gaussian distribution.

Figures 4 and 5 show how voltage and temperature effect the mean and standard deviation of the BMX055 z-axis acceleration measurement. Temperature has a greater effect on mean and standard deviation than voltage. Both voltage and temperature have a greater effect upon the standard deviation than the mean. A random number generator based on this phenomenon should control both the temperature and the supply voltage of the sensor.

Fig. 4: The mean of BMX055 z-axis distributions increases with increasing temperature and supply power. Normalised to the maximum observed mean.
Fig. 5: The standard deviation of BMX055 z-axis distributions decreases with increasing temperature and supply power. Normalised to the maximum observed standard deviation.

Figure 6 shows the effect of increasing the bin size on the divergence between a distribution of 100,000 values and its ideal fitted distribution. This shows that increased quantization increases the difference between a distribution and its fitted Gaussian.

Fig. 6: The effect of increasing discretisation on the KL divergence whilst keeping the total number of samples constant at 100,000. Increased quantisation increases KL divergence.

V Gaussian To Gaussian Transform

A univariate Gaussian can be transformed to any other univariate Gaussian with one multiplication and one addition. This is significantly less computation than the accept reject method which requires at least 10 operations per each accept reject test: an exponential, square, square root, subtraction, comparison and five divisions / multiplications. The accept reject method may furthermore need to repeat this set of operations numerous times for each random variate. Two uniform random numbers must be generated for one accept reject method sample. Rapid addition and multiplication can be achieved using fast adders and multipliers implemented on an FPGA, Figure 7 shows how this could be achieved. The CPU requests a distribution by specifying parameters to the transform circuitry which proceeds to transform the input Gaussian to fit the requested output distribution. The transform circuitry then stores the values in a small high-speed cache that the CPU can read from. The asynchronous transform circuitry on the FPGA will be free to execute transformations at a rate not bound by the global clock frequency. Offloading the transformation to dedicated hardware leaves the processor free to execute other instructions which will improve performance.

Fig. 7: Gaussian to Gaussian transform implementation. The CPU can request distributions with certain parameters and they are presented to it in a small fast cache.

Vi Conclusion

Sensors are a feasible source of non-uniform random variates at a higher sample rate and with greater efficiency than all of the state-of-the-art methods used to generate them. The mean and standard deviation of the noise produced by the z-axis of a commercial accelerometer depend upon the temperature of the environment and the supply voltage. The parameters of the distribution drift with temperature and supply voltage so both must be kept constant to avoid this. Quantizing a Gaussian distribution increases the KL divergence between it and a fitted Gaussian. This type of non-uniform random number generator can reduce the error of Monte Carlo integration of a univariate Gaussian by a factor of whilst doubling the speed.


  • [1] E. Aamir, Z. K. Nagy, C. D. Rielly, T. Kleinert, and B. Judat (2009-09-16)

    Combined quadrature method of moments and method of characteristics approach for efficient solution of population balance models for dynamic modeling and crystal size distribution control of crystallization processes

    Industrial & Engineering Chemistry Research 48 (18), pp. 8575–8584. Cited by: §I-C.
  • [2] M. Chang (2010) Monte carlo simulation for the pharmaceutical industry: concepts, algorithms, and case studies. CRC Press. External Links: ISBN 9781439835920 Cited by: §I-C.
  • [3] R. L. Cook (1986-01) Stochastic sampling in computer graphics. ACM Trans. Graph. 5 (1), pp. 51–72. External Links: ISSN 0730-0301 Cited by: §I-C.
  • [4] Y. Dai, W. Jiang, and G. Wang (2016-08) Building bayesian inference graphs for healthcare statistic evidence. In 2016 45th International Conference on Parallel Processing Workshops (ICPPW), Vol. , pp. 415–420. External Links: ISSN 2332-5690 Cited by: §I-C.
  • [5] J. Danger, A. Ghazel, E. Boutillon, and H. Laamari (2000-12) Efficient fpga implementation of gaussian noise generator for communication channel emulation. In ICECS 2000. 7th IEEE International Conference on Electronics, Circuits and Systems (Cat. No.00EX445), Vol. 1, pp. 366–369 vol.1. External Links: ISSN Cited by: §I-C.
  • [6] L. Devroye (1986) Non-Uniform Random Variate Generation. Springer-Verlag, McGill University Montreal H3A 2K6 Canada. Cited by: §I-B.
  • [7] KLDIV. Nima Razavi. Note: [Online]. Available: https://uk.mathworks.com/matlabcentral/fileexchange/13089-kldiv?s_tid=FX_rc1_behav, Accessed: 27/03/2019 Cited by: §IV.
  • [8] D. P. Kroese, T. Brereton, T. Taimre, and Z. I. Botev (2014) Why the monte carlo method is so important today. Wiley Interdisciplinary Reviews: Computational Statistics 6 (6), pp. 386–392. Cited by: §I-C.
  • [9] B. Lambert

    A student’s guide to bayesian statistics

    SAGE. Cited by: §I-C, §I-C.
  • [10] V. V. Marcel, I. J.B.F. Adan, and S. A.E. Resing-Sassen (2006) Stochastic modeling of manufacturing systems. Springer. External Links: ISBN 978-3-540-29057-5 Cited by: §I-C.
  • [11] M. Mitzenmacher and E. Upfal (2005) Probability and computing: randomized algorithms and probabilistic analysis. Cambridge University Press, New York, NY, USA. External Links: ISBN 0521835402 Cited by: §I-C.
  • [12] J. R Oaks, K. A Cobb, V. N Minin, and A. D Leaché (2019-01) Marginal Likelihoods in Phylogenetics: A Review of Methods and Applications. Systematic Biology 68 (5), pp. 681–697. Cited by: §I-C.
  • [13] Randomness requirements for security. Network Working Group. Cited by: §I-A.
  • [14] M. N. O. Sadiku and M. Ilyas (1995) Simulation of local area networks. CRC Press, Inc., Boca Raton, FL, USA. External Links: ISBN 0849324734 Cited by: §I-C.
  • [15] D. B. Thomas and W. Luk (2007-07) Non-uniform random number generation through piece-wise linear approximations. IET Computers Digital Techniques 1 (4), pp. 312–321. External Links: ISSN 1751-8601 Cited by: §I.
  • [16] D. B. Thomas and W. Luk (2008-07)

    Resource efficient generators for the floating-point uniform and exponential distributions

    In 2008 International Conference on Application-Specific Systems, Architectures and Processors, Vol. , pp. 102–107. External Links: ISSN 1063-6862 Cited by: §I-C.
  • [17] D.B. Thomas and W. Luk (2006-04) Efficient hardware generation of random variates with arbitrary distributions. In 14th Annual IEEE Symposium on Field-Programmable Custom Computing Machines, Vol. , pp. 57–66. External Links: ISSN Cited by: §I-B, §I.
  • [18] D. B. Thomas, L. Howes, and W. Luk (2009) A comparison of cpus, gpus, fpgas, and massively parallel processor arrays for random number generation. In Proceedings of the ACM/SIGDA International Symposium on Field Programmable Gate Arrays, FPGA ’09, New York, NY, USA, pp. 63–72. External Links: ISBN 978-1-60558-410-2 Cited by: §I-A, TABLE I, §I.
  • [19] S. Wang, A. R. Lebeck, and C. Dwyer (2015-Sep.) Nanoscale resonance energy transfer-based devices for probabilistic computing. IEEE Micro 35 (5), pp. 72–84. External Links: ISSN 0272-1732 Cited by: §I.
  • [20] G. L. Zhang, P. H. W. Leong, C. H. Ho, K. H. Tsoi, C. C. C. Cheung, D. Lee, R. C. C. Cheung, and W. Luk (2005-12) Reconfigurable acceleration for monte carlo based financial simulation. In Proceedings. 2005 IEEE International Conference on Field-Programmable Technology, 2005., Vol. , pp. 215–222. External Links: ISSN Cited by: §I-C.