Pitfalls in pseudo-random number sampling at scale with Apache Spark

With kind contributions from Gregory Piatetsky-Shapiro, KDnuggets

In many data science applications and in academic research, techniques involving Bayesian Inference is now used commonly. One of the basic operation in Bayesian Inference techniques is drawing instances from given statistical distribution. This of course well known pseudo-random number sampling. Most commonly used methods first generates uniform random number and mapping that into distribution of interest via cumulative function (CDF) of it, i.e., Box-Mueller transform.

Large scale simulation are now possible, due to highly stable computational frameworks that can scale well. One of the unique framework is Apache Spark due to its distributed data structure supporting fault tolerance, called Resilient Distributed Data (RDD). Here is a simple way to generate one million Gaussian Random numbers and generating an RDD:

1 2 3 4 5 6

// Generate 1 million Gaussian random numbers import util.Random Random.setSeed(4242) val ngauss = (1 to 1e6.toInt).map(x => Random.nextGaussian) val ngauss_rdd = sc.parallelize(ngauss) ngauss_rdd.count // 1 million

One unrealistic part of the above code example is that you may want to generate huge number of samples that won't fit in single memory, ngauss variable above.  Luckily, there are set of library functions one can use to generate random data as an RDD from mllib, see randomRDD. But for the remainder of this post, we will use our home made random RDD.

Figure: Scaling of execution time with increasing size,
with or without re-partitioning.

Concept of Partitions

As RDDs are distributed data structures, the concept of partition comes into play (link)..  So, you need to be careful of the size of partitions in RDDs. Recently I posed a question about this in Apache Spark mailing list (link)(gist). If you reduce the data size, take good care that your partition size reflects this, so to speak avoiding huge performance reduction. Unfortunately, Spark does not provide an automated  out of box solution optimising partition size. The actual data items that might reduce during your analysis pipeline. A reduced RDD will inherit partition size of its parent and this may be a limiting issue.

As you might have already guessed, RDDs are great tool in doing large scale analysis but they won't provide you a free lunch. Let's do a small experiment.

Hands on Experiment

Going back to our original problem of using Spark in Bayesian inference algorithms, it is common to operate on samples via certain procedure. And those procedures, let's say an operator, highly likely that it will reduce the number of elements in the sample. One example would be applying a cut-off or a constrained in the CDF, which essential the definition of it, probability of random variable $x > x_{0}$.  As seen in Figure, we have generated random RDDs up to 10 million numbers and measure the wall-clock time of `count` operation,  which occurs after a filter operation that reduces the number of items considerably. See codes in the Appendix. As a result, in Figure, we have identified 3 different regions, depending on data size, 

1.  Small Data: Re-partitioning does not play a role in the performance.
2.  Medium Size: Re-partitioning gives up to order of magnitude better performance.
3.  Big Data: Re-partitioning gives a constant performance improvement, up to 3 times better, and  the improvement is drifting, meaning it will be more significant larger the data size. 

Conclusion

Spark provides a superb API to develop high quality Data Science solutions. However, programming with Spark and designing algorithms requires optimisation of different aspects of the RDD workflow. Here, we only demonstrate only dramatic effect of re-partitioning after a simple operation in the walk clock time. Hence, it is advised to have a benchmark identifying under which circumstances your data pipeline produce different wall clock behaviour before going into production.

Appendix

Entire code base can be cloned from github  (here).

Spark Benchmark

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

/* Generating Gaussian numbers : Performance with repartition RDD partition size retains from the parent (c) 2017 GPLv3 Author: Mehmet Suzen (suzen at acm dot org) Run this script on spark-shell spark-shell --executor-cores 4 spark-shell> :load gaussian_random_filter.scala */ import util.Random // normal import breeze.linalg._ // Linear Algebra objects and csvwrite import java.io.File // File io /* Generate gaussian random numbers manually without mllib Benchmark repartition */ // Random numbers to generate val Ns = Array(1e3, 1e4, 5e4, 1e5, 5e5, 1e6, 2e6, 4e6, 8e6, 1e7) val benchMat = DenseMatrix.zeros[Double](10,3) Random.setSeed(4242) for(i <- 0 to Ns.size-1) { println("running for " + Ns(i)) // Generate random RDD size Ns var ngauss = (1 to Ns(i).toInt).map(x=>Random.nextGaussian) var ngauss_rdd = sc.parallelize(ngauss) var ngauss_rdd2 = ngauss_rdd.filter(x=>x > 4.0) // An operation without repartition var t0 = System.nanoTime() var c1 = ngauss_rdd2.count var t1 = System.nanoTime() var e1 = (t1 - t0)/1e9 // seconds // An operation with repartition var ngauss_rdd3 = ngauss_rdd2.repartition(1) t0 = System.nanoTime() var c2 = ngauss_rdd3.count t1 = System.nanoTime() var e2 = (t1 - t0)/1e9 benchMat(i,::) := DenseVector[Double](Ns(i), e1, e2).t } /* Record the benchmark results */ csvwrite(new File("bench.csv"), benchMat)

Plotting code

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

# # Plot the benchmark from Spark Gaussian Random numbers # # (c) 2017 # GPLv3 # library(grid) library(gridExtra) library(ggplot2) library(reshape) bench_df <- read.csv2("bench.csv", header=FALSE, sep=",") colnames(bench_df) <- c("N", "no", "yes") bench_df2 <- reshape2:::melt(bench_df, measure.vars=c("no","yes")) colnames(bench_df2) <- c("N", "repartion", "time") bench_df2$N <- as.numeric(as.vector(bench_df2$N)) bench_df2$time <- as.numeric(as.vector(bench_df2$time)) gt <- theme( panel.background = element_blank(), axis.text.x = element_text(face="bold", color="#000000", size=11), axis.text.y = element_text(face="bold", color="#000000", size=11), axis.title.x = element_text(face="bold", color="#000000", size=11), axis.title.y = element_text(face="bold", color="#000000", size=11) ) p1 <- ggplot(bench_df2, aes(x=N, y=time, colour=repartion)) + geom_smooth(formula="y~x", span=0.3) + xlab("Number of random draws") + ylab("Wall Clock (Seconds)") + ggtitle("Effect of repartition in count: Gaussian Random Numbers") + gt grid.newpage() footnote <- "(c) 2017, Mehmet Suzen : http://memosisland.blogspot.de/" g <- arrangeGrob(p1, bottom = textGrob(footnote, x = 0, hjust = -0.1, vjust=0.1, gp = gpar(fontface = "italic", fontsize = 12))) png(file="spark_repartition_random.png") grid.draw(g) dev.off()

S-shaped data: Smoothing with quasibinomial distribution

Figure 1: Synthetic data and fitted curves.

S-shaped distributed data can be found in many applications. Such data can be approximated with logistic distribution function [1].  Cumulative distribution function of logistic distribution function is a logistic function, i.e., logit.

To demonstrate this, in this short example, after generating a synthetic data, we will fit quasibinomial regression model to different observations.

ggplot [2], an implementation of grammar of graphics [3], provides capability to apply regression or customised smoothing onto a raw data during plotting.

Generating Synthetic Data

Let generate set of $n$ observation over time $t$, denoted, $X_{1}, X_{2}, ..., X_{n}$ for $k$ observation $X=(x_{1}, x_{2}, ..., x_{k})$. We will use cumulative function for logistic distribution [4],
$$F(x;\mu,s) = \frac{1}{2} + \frac{1}{2} \tanh((x-\mu)/2s)$$, adding some random noise to make
it realistic.

Let's say there are $k=6$ observations with the following parameter sets, $\mu = \{9,2,3,5,7,5\}$  and $s=\{2,2,4,3,4,2\}$, we will utilise `mapply` [5] in generating a syntetic data frame.

generate_logit_cdf <- function(mu, s, 
                               sigma_y=0.1, 
                               x=seq(-5,20,0.1)) {
  x_ms <- (x-mu)/s
  y    <- 0.5 + 0.5 * tanh(x_ms)  
  y    <- abs(y + rnorm(length(x), 0, sigma_y))
  ix   <- which(y>=1.0)
  if(length(ix)>=1) { 
    y[ix] <- 1.0
  }
  return(y)
}
set.seed(424242)
x      <- seq(-5,20,0.025) # 1001 observation
mu_vec <- c(1,2,3,5,7,8)   # 6 variables
s_vec  <- c(2,2,4,3,4,2)
# Syntetic variables
observations_df<- mapply(generate_logit_cdf, 
                              mu_vec, 
                              s_vec, 
                              MoreArgs = list(x=x))
# Give them names
colnames(observations_df) <- c("Var1", "Var2", "Var3", "Var4", "Var5", "Var6")
head(observations_df)  

Smoothing of observations

Using the syntetic data we have generated, `observations_df`,
we can noq use `ggplot` and `quasibinomial` `glm` to visualise
and smooth the variables.

library(ggplot2)
 library(reshape2)
 df_all <- reshape2:::melt(observations_df)
 colnames(df_all) <- c("x", "observation", "y")
 df_all$observation <- as.factor(df_all$observation)
 p1<-ggplot(df_all, aes(x=x, y=y, colour=observation)) + geom_point() + 
        scale_color_brewer(palette = "Reds") +
        theme(
              panel.background = element_blank(), 
              axis.text.x      = element_text(face="bold", color="#000000", size=11),
              axis.text.y      = element_text(face="bold", color="#000000", size=11),
              axis.title.x     = element_text(face="bold", color="#000000", size=11),
              axis.title.y     = element_text(face="bold", color="#000000", size=11)
#              legend.position = "none"
              )
 l1<-ggplot(df_all, aes(x=x, y=y, colour=observation)) +
        geom_point(size=3) + scale_color_brewer(palette = "Reds") + 
        scale_color_brewer(palette = "Reds") + 
        #geom_smooth(method="loess", se = FALSE, size=1.5) +
        geom_smooth(aes(group=observation),method="glm", family=quasibinomial(), formula="y~x",
                       se = FALSE, size=1.5) +
        xlab("x") + 
        ylab("y") +
        #scale_y_continuous(breaks=seq(0.0,1,0.1)) +
        #scale_x_continuous(breaks=seq(0.0,230,20)) +
        #ggtitle("")  + 
        theme(
              panel.background = element_blank(), 
              axis.text.x      = element_text(face="bold", color="#000000", size=11),
              axis.text.y      = element_text(face="bold", color="#000000", size=11),
              axis.title.x     = element_text(face="bold", color="#000000", size=11),
              axis.title.y     = element_text(face="bold", color="#000000", size=11)
              )
 library(gridExtra)
 gridExtra:::grid.arrange(p1,l1)

References
[1] https://en.wikipedia.org/wiki/Logistic_distribution#Applications
[2] http://www.ggplot.org
[3] The Grammar of Graphics, L. Wilkinson, http://www.amzn.com/038724544
[4] http://en.wikipedia.org/wiki/Logistic_distribution#Cumulative_distribution_function.
[5] https://stat.ethz.ch/R-manual/R-devel/library/base/html/mapply.html

Notes:
* Quasibinomial distribution is as a term used in R's GLM implementation context, see here.

Euclid Algorithm for Set of Integers: 'Reduce' vs. trees in R

The Euclid Algorithm provides a solution to the greatest common divisor (GCD) of two natural numbers $x_{1}$ and $x_{-2}$, denoted by $GCD(x_{1}, x_{2})$. This will produce the largest integer that divides $x_{1}$ and $x_{2}$. Solution is proposed by Euclid of Ancient Greece. This can be formulated as a recurrence relation:
$$x_{k} = x_{k-1}  modulo  x_{k-2}.$$
'modulo' binary operation returns the remainder from a given division. Stopping criterion for the recurrence relation is reached when $x_{k-2}=0$ and the result of GCD will be the current value of $x_{k-1}$. This process can be visualised as successive divisions. Let's implement this in R in a naive way.

# Naive Euclid algorithm by msuzen
gcd <- function(a, b)
{
  rk_1 <- a;
  rk_2 <- b;
  # Recurrence Formula:  r_k =  r_k-1 modulo r_k-2
  # Increment k until r_k-2 == 0 
  while(rk_2 != 0) {
    rk      <- rk_1%%rk_2; # remainder
    rk_1    <- rk_2;       # proceed in recurrence
    rk_2    <- rk;
  }
  return(rk_1)
}

This is a straight forward task. Let's make the problem little more generic. What happends if we would like to know GCD of $n$ natural numbers, $x_{1},..., x_{n}$? Than, a solution is to apply GCD operation pairwise, for example if $n=3$:
$$GCD(x_{1}, GCD(x_{2}, x_{3})) = GCD(GCD(x_{1}, x_{2}), x_{3})$$
How can we implement this for a vector of non-negative integers?

Tree Approach

The simplest way to reach GCD of $n$ numbers is probably thinking of this process as a binary tree, formed by pairing elements of set of integers as we obtain GCDs. It is relatively easy to implement this because ordering of pairs is not important. We can start from the beginning and pair up as we obtain the results. Here is the naive implementation.

1 2 3 4 5 6 7 8

gcdN <- function(X) { n <- length(X) gcdEnd <- X[1] for(i in 2:n) { gcdEnd <- gcd(gcdEnd, X[i]) } return(gcdEnd) }

'Reduce' operation 

So called tree approach we have given above is actually noting but a Reduce operation in the context of MapReduce. The function gcdN can be replaced with a single line.

Reduce("gcd", X)

This example looks trivial but having Reduce and friends in our programming toolbox makes our life a little easier. Noting so novel here! But it could saves us time not to implement tree like algorithm.

Metaprogramming in R with an example: Beating lazy evaluation

Functional languages allows us to treat functions as types. This brings us a distinct advantage of being able to write a code that generates further code, this practise is generally known as metaprogramming. As a functional language R project provides tools to perform well structured code generation. In this post, I will present a simple example that generates functions on the fly based on different parametrisation in the function body. Consider the following simple function taking a vector as an argument and returning the number of element that are higher than a given threshold. 

1 2 3 4

myFun <- function(vec) { numElements <- length(which(vec > threshold)) numElements }

If somehow we need to have a different threshold value within the body, for a moment accept that it is a requirement rather than proposing to have an other argument in the function definition. Instead of rewriting the function by hand we will write a function that generates all these functions in our work space. Problematic bit of this exercise will be to beat lazy evalution.  Here is the function that produces losts of myFun type functions:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

genMyFuns <- function(thresholds) { ll <- length(thresholds) print("Generating functions:") for(i in 1:ll) { fName <- paste("myFun.", i, sep="") print(fName) assign(fName, eval( substitute( function(vec) { numElements <- length(which(vec > tt)); numElements; }, list(tt=thresholds[i]) ) ), envir=parent.frame() ) } }

Let's shortly analyse  this function. If we don't use substitute explicitly there, due to lazy evalution our value for the threshold will not be assigned at the loop value but the  last value of thresholds[i]. Here is one numeric example on the R CLI session:

1 2 3 4 5 6 7 8 9 10 11 12

> genMyFuns(c(7, 9, 10)) [1] "Generating functions:" [1] "myFun.1" [1] "myFun.2" [1] "myFun.3" > myFun.1(1:20) [1] 13 > myFun.2(1:20) [1] 11 > myFun.3(1:20) [1] 10 >

To be able to generate code is very powerful tool. However, a caution should be taken in practicing code generation in a large project. This may bring more problems in debugging. Every powerful method comes with a hidden cost.

Practicing static typing in R: Prime directive on trusting our functions with object oriented programming

The creator of S language which R is derived from John Chambers said in one of his books  Software for data analysis programming with R: 

...This places an obligation on all creators of software to program in such a
way that the computations can be understood and trusted. This obligation
I label the Prime Directive.

He was referring to prime directive from Star Trek. One of the practice in this direction is to have a proper checks in place for the types we use. We can trust that if we pass for example a wrong type to our function, it will fail gracefully. So a type system of a programming language is quite important in mission critical numerical computations. Since R language is weakly typed language, or dynamically typed similar to Perl, Python or Matlab/Octave, most of R users omit to place type checks in their functions if not rarely. For example take the following function that takes arguments of a matrix, a vector and a function name. It applies the named function to each columns of the matrix listed in the given vector. Assuming named function is returning a single number our function will return a vector of numbers.

myMatrixOperation  <-  function(A, v, fName) {
  sliceA <-  A[, v];
  apply(sliceA, 2, fName);
}

One obvious way to put if statements for each argument in our function. So, function may look like:

myMatrixOperation <- function(A, v, fName) {
  if(!is.matrix(A)) {
   stop("A is not a matrix");
  }
  if(!is.vector(v)) {
   stop("v is not a vector");
  }
  if(!is.funcion(fName)) {
   stop("fName is not a function");
  }
  sliceA <- A[, v];
  apply(sliceA, 2, fName);
}

The problem with this approach appears to be the fact that it is too verbose and if we have a repeating pattern of arguments in many functions and many arguments, we would copy and paste code many times. It would not only look ugly but wastes our time. Luckily there is a mechanism to address this: S4-class system. Let's define an S4 class for our set of arguments, following an example instantiation.

setClass("mySlice", representation(A="matrix", v="vector", fName="function"))

myS <- new("mySlice",, A=matrix(rnorm(9),3,3),v=c(1,2), fName=mean)
str(myS)
Formal class 'mySlice' [package ".GlobalEnv"] with 3 slots
  ..@ A    : num [1:3, 1:3] 0.356 -0.34 -0.642 -0.466 2.915 ...
  ..@ v    : num [1:2] 1 2
  ..@ fName:function (x, ...)

Now if we re-write the function that uses our S4 class with type checking only to passing object once.

is.mySlice <- function(obj) {
  l <- FALSE 
  if(class(obj)[1] == "mySlice") { l <- TRUE } 
  l 
} 

myMatrixOperation <- function(mySliceObject) { 
  if(!is.mySlice(mySliceObject)) { 
    stop("argument is not class of mySlice") 
  }  
  sliceA <- mySliceObject@A[, mySliceObject@v]; 
  apply(sliceA, 2, mySliceObject@fName); 
} 

This simple example demonstrates how we can introduce a good organization to our R codes, that obeys the prime directive. Further more modern approach to object orientation is introduced by John Chambers called Reference classes. If you practice this kind of approach in your R codes than I can only say; Live long and prosper.

Function closures and abstraction in MATLAB/Octave

Recently a short idea on how function closures can be build in R is shown [link]. It is also discussed in detail by Darren Wilkinson's blog [link]. Here we will mimic the approach presented there. Recall that function closures are basically a reference to a function with the attachment of its lexical environment. Note that, it is not a function pointer while it also return the values of its arguments if they are defined in the given scope. This sort of language construction can be very handy when some of the inner details need to be parametrized, so to speak to prevent re-writing the function over and over again.

A trivial example of adding a number to a given number can be formulated.  Let's say x is a number and we need to add 4 to x.

function add4(x)
  x+4
end

So far so good, but do you really want to re-write an other function to add 7. Not really! So instead we can generate a function that generates addition function for 7 (Sounds like an onion.). This is quite standard in functional languages and noting new actually!

>> addx = @(x) @(y) x + y ;
>> add4=addx(4);
>> add7=addx(7);
>> add4(5)
9 
>> add7(5)
12

So function add4 and add7 are generated by the function addx. We just utilized function handles. This might be quite trivial. Let's try the technique in a little more realistic setting. Imagine we need to take a median of each vector, which has different sizes stored in a structure.  Let's generate one

>> rng(0.42,'v4'); % reproduce numbers
>> myStr.vec1 =  rand(11,1);
>> myStr.vec2 =  rand(5,1);
>> myStr.vec3 =  rand(20,1);
>> myStr.vec4 =  rand(8,1);

So a function that takes a structure like this and returns the median value of the named member

>> medianVec= @(str, memName) median(str.(memName));

For example for vec1

>> medianVec(myStr, 'vec1');
 0.5194

If we need this procedure for vec1 repetitively, we don't need to re-write the second argument all the time. So

>> medianVecN = @(memName)  @(str) median(str.(memName));
>> medianVec1 = medianVecN('vec1');
>> medianVec1(myStr); % Vala!
0.5194 

This might still look trivial but hopefully it would give you an idea of the technique.