\documentclass{beamer}  
\usepackage{graphicx} 
\usepackage{url}
\usepackage{fancyvrb}
\usepackage{relsize}
\usepackage{listings}

\AtBeginSection[]
{
\begin{frame}<beamer>{Table of Contents}
\tableofcontents[currentsection,currentsubsection, 
    hideothersubsections, 
    sectionstyle=show/shaded,
]
\end{frame}
}


% \lstset{language=R}

\mode<presentation>
{ \usetheme{Hannover} }

\title{Data Science Applications of GPUs in the R Language}
\author{Norm Matloff \\
University of California at Davis 
\ \\
\ \\
GTC 2016
}

\date{April 7, 2016}

\begin{document} 

% \lstset{language=R}
\lstset{language=R,basicstyle=\smaller,commentstyle=\rmfamily\smaller,
   showstringspaces=false,
   xleftmargin=4ex,literate={<-}{{$\leftarrow$}}1 {~}{{$\sim$}}1}

\begin{frame}
\titlepage

% URL for these slides (repeated on final slide):
These slides at
\url{http://heather.cs.ucdavis.edu/GTC.pdf}

\end{frame}

\begin{frame}
\frametitle{Why R?}
\pause

\begin{itemize}

\item The {\it lingua franca} for the data science
community.  (R-Python-Julia battle looming?)

\item Statistically Correct: Written by statisticians, for statisticians.

\item 8,000 CRAN packages!

\item Excellent graphics capabilities, including Shiny (easily build your own 
interactive tool).

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{R $\rightarrow$ GPU Link Pros and Cons}
\pause

On the plus side:

\begin{itemize}

\item Speed: R is an interpreted language.  (Nick Ulle and Duncan Temple
Lang working on LLVM compiler.)

\item R is often used on large and/or complex data sets, thus requiring
large amounts of computation.

\item Much of R computation involves matrices or other operations
well-suited to GPUs.

\end{itemize}

On the other hand:

\begin{itemize}

\item Big Data implies need for multiple kernel calls, and much
host/device traffic.

\item Ditto for R's many iterative algorithms.

\item Many of the matrix ops are not embarrassingly parallel.

\item Unpacking and repacking into R object structure.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Disclaimers}\pause

\begin{itemize}

\item Talk is meant to be aimed at NVIDIA but otherwise generic, not
focusing on the latest/greatest model.
\pause

\item Our running example, NMF, has the goal of illustrating issues and
methods concerning the R/GPU interface. It is not claimed to produce the
fastest possible computation.  (See talk by Wei Tan in this session.)

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Running Example: Nonnegative Matrix Factorization (NMF)}\pause

\begin{itemize}

\item  Have matrix $A \geq 0$, rank $r$.

\item Want to find matrices $W \geq 0$ and $H \geq 0$ of rank $s \ll r$
with

$$
A \approx W H
$$

\item Columns of $W$ form a ``pseudo-basis'' for columns of A:
$A_{.j}$ is approximately a linear combination of the columns of $W$,
with coordinates in $H_{.j}$.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Applications of NMF}\pause

\begin{itemize}

\item Image compression.
\pause

\item Image classification. \pause Each column of $A$ is one image.
\pause To classify new image, find coordinates $u$ w.r.t.\ $W$, then
find nearest neighbor(s) of u in $H$.
\pause

\item Text classification.  Each column of $A$ is one document, with
counts of words of interest. Similar to image classification.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Example of R Calling C/C++}\pause

\begin{itemize}

\item Compare R's {\bf NMF} package to E. Battenberg's {\bf NMF-CUDA},
on a 3430 $\times$ 512 $A$:
\pause

\item R, $s = 10$:  649.843 sec

\item GPU, $s = 30$:  0.986 sec

\item GPU solved a much bigger problem in much less time

\item Even though the R pkg is in C++, not R.
\pause

\item Solution: Call {\bf NMF-CUDA}'s {\bf update\_div()} from R.
\pause BUT HOW?

\item R's {\bf Rcpp} package makes interfacing R to C/C++ very
convenient and efficient.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{General R/GPU Tools}
\pause

What's out there now for R/GPU:

\begin{itemize}

\item {\bf gputools} 

(Buckner {\it et al.}) The oldest major package.  Matrix multiply; matrix
of distances between rows; linear model fit; QR decomposition;
correlation matrix; hierarchical clustering.

\item {\bf HiPLAR}

(Montana {\it et al}.)
R wrapper for {\bf MAGMA} and {\bf PLASMA}.  Linear algebra routines,
e.g. Cholesky.

\item {\bf rpud}

(Yau.) Similar to {\bf gputools}, but has SVM.

\item {\bf Rth}

(Matloff.)  R interfaces to some various algorithms coded in Thrust.
Matrix of distances between rows; histogram; column sums; Kendall's Tau;
contingency table.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Current Tools (cont'd.)}
\pause

\begin{itemize}

\item {\bf gmatrix}

(Morris.)
Matrix multiply, matrix subsetting, Kronecker product, row/col sums, 
Hamiltonian MCMC, Cholesky.  

\item {\bf RCUDA}

(Baines and Temple Lang, currently not under active development.)
Enables calling GPU kernels directly from R.  (Kernels still written in
CUDA.) 

\item {\bf rgpu}

(Kempenaar, no longer under active development.) ``Compiles'' simple
expressions to GPU.

\item various OpenCL interfaces

{\bf ROpenCL}, {\bf gpuR}.  Similar to {\bf RCUDA}, but via OpenCL
interface.

\end{itemize}

\end{frame}

\begin{frame}[fragile]
\frametitle{Example:  Linear Regression Via gputools}\pause

\begin{lstlisting}
> test <- function(n,p) {
   x <- matrix(runif(n*p),nrow=n)
   regvals <- x %*% rep(1.0,p)
   y <- regvals + 0.2*runif(n)
   xy <- cbind(x,y)
   print("gputools method")
   print(system.time(gpuLm.fit(x,y)))
   print("ordinary method")
   print(system.time(lm.fit(x,y)))
}
> test(100000,1500)
[1] "gputools method"
   user  system elapsed
  6.280   2.878  17.902
[1] "ordinary method"
   user  system elapsed
142.282   0.669 142.912
\end{lstlisting}


\end{frame}

\begin{frame}
\frametitle{Key Issue: Keeping Objects on the Device}
\pause

\begin{itemize}

\item Some packages, notably {\bf gputools}, do not take arguments on
the device.

\item So, cannot store intermediate results on the device, thus requiring
needless copying.

\item Some packages remedy this, e.g. {\bf gmatrix}.
 
\end{itemize}

\end{frame}

\begin{frame}[fragile]
\frametitle{Example}
\pause

\begin{lstlisting}
library(gputools)
library(gmatrix)
n <- 5000
z <- matrix(runif(n^2),nrow=n)
# plain R:
system.time(z %*% z %*% z)
#    user  system elapsed 
# 138.757   0.322 139.081 
system.time(gpuMatMult(gpuMatMult(z,z),z))
#    user  system elapsed 
#   6.607   1.170  10.059 
zm <- gmatrix(z,nrow=n,ncol=n) # zm2, zm3 not shown
system.time({gmm(zm,zm,zm2); gmm(zm,zm2,zm3)})
#    user  system elapsed 
#   6.258   1.031   7.285 
\end{lstlisting}

\end{frame}

\begin{frame}
\frametitle{Rth Example --- Kendall's Tau}\pause

A kind of correlation measure, defined to be the proportion of {\it
concordant pairs}:

$(X_i,Y_i) \textrm{ and } (X_j,Y_j)$ are concordant if
$sign(X_i-X_j) \cdot sign(Y_i-Y_j) > 0$

\end{frame}

\begin{frame}[fragile]
\frametitle{Kendall's Tau (cont'd.)}\pause 

R wrapper to Thrust call:

\begin{lstlisting}
rthkendall <- function(x,y) {
   dyn.load("rthkendall.so")
   n <- length(x)
   tmp <-
      .C("rthkendall",as.single(x),as.single(y),
      as.integer(n),tmpres=single(1),DUP=dupval)
   return(tmp$tmpres)
}
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]
\frametitle{Kendall's Tau (cont'd)}\pause

\begin{lstlisting}
void rthkendall(float *x, float *y, 
   int *nptr, float *tauptr)
{  int n = *nptr;
   thrust::counting_iterator<int> seqa(0);
   thrust::counting_iterator<int> seqb = seqa + n-1;
   // dx, dy, tmp declarations not shown
   thrust::transform(seqa,seqb,tmp.begin(),
      calcgti(dx,dy,n));
   int totcount = 
      thrust::reduce(tmp.begin(),tmp.end());
   float npairs = n * (n-1) / 2;
   *tauptr = (totcount - (npairs-totcount)) / npairs;
}
\end{lstlisting}

\end{frame}

\begin{frame}[fragile]
\frametitle{Kendall's Tau (cont'd)}\pause

\begin{lstlisting}
struct calcgti {  // handle 1 i, all j > i
   // more declarations not shown
   calcgti(floublevec _dx,floublevec _dy,int _n) :
      dx(_dx),
      dy(_dy),
      n(_n)
      {  wdx = thrust::raw_pointer_cast(&dx[0]);   
         wdy = thrust::raw_pointer_cast(&dy[0]);  }
   __device__ int operator()(int i)
   {  flouble xi = wdx[i], yi = wdy[i];
      int j,count=0;
      for (j = i+1; j < n; j++)
         count += 
            ( (xi - wdx[j]) * (yi - wdy[j]) > 0);
      return count;
   }
};
\end{lstlisting}

\end{frame}

\begin{frame}
\frametitle{Example: NMF Again}\pause

\begin{itemize}

\item The R {\bf NMF} package, and {\bf NMF-CUDA} use {\it
multiplicative update} methods.

\item For instance, for Frobenius norm,

$$
W \leftarrow W \circ 
\frac
{AH'}
{WHH'}
$$

% $$
% H_{kj} \leftarrow H_{kj}
% \frac
% { \sum_{i} W_{ik} A_{ij} / (WH)_{ij}}
% {\sum_{i} W_{ik}}
% $$

and similarly for $H$.
\pause

\item Another possibility is to use the {\it alternating least squares}
method:

   \begin{itemize}

   \item In odd-numbered iterations, regress each col. of $A$ against
   cols. of $W$, yielding the columns of $H$.  Mult. update even better
   suited to GPUs.

   \item In even-numbered iterations, reverse the roles of $W$ and $H$
   (and now with rows).

   \end{itemize}

\item As seen earlier, least-squares estimation can be done fairly well
on GPUs.
\end{itemize}

\end{frame}

\begin{frame}[fragile]
\frametitle{RCUDA Example: Normal Density}\pause

Basic goal: Call CUDA kernels from R without burdening the R programmer
with details of configuring grids, allocating device memory, copying
between host and device, etc.

Kernel:

\begin{lstlisting}
extern "C"
__global__ void 
  dnorm_kernel(float *vals,int n,float mu,float sig)
{
  int myblock = blockIdx.x + blockIdx.y * gridDim.x;
  int blocksize = 
     blockDim.x * blockDim.y * blockDim.z;
  int subthread = 
     threadIdx.z*(blockDim.x * blockDim.y) + 
     threadIdx.y*blockDim.x + threadIdx.x;
  int idx = myblock * blocksize + subthread
  float std = (vals[idx] - mu)/sig;
  float e = exp( - 0.5 * std * std);
  vals[idx] = e / ( sig * sqrt(2 * 3.14159));
}

\end{lstlisting}

\end{frame}

\begin{frame}[fragile]
\frametitle{RCUDA (cont'd.)}\pause
\begin{lstlisting}
n = 1e6
mean = 2.3
sd =  2.1
x = rnorm(n, mean, sd)
# eval density at all pts in x
m = loadModule("dnorm.ptx")
k = m$dnorm_kernel
ans = .cuda(k,x,n,mean,sd, 
   gridDim = c(62, 32),blockDim = 512)
\end{lstlisting}

\end{frame}

\begin{frame}
\frametitle{Helpful Utilities}
\pause

\begin{itemize}

\item {\bf Rcpp}

   \begin{itemize}

   \item Greatly facilitates calling C/C++ from R.

   \item Base R offers functions {\bf .C()} and {\bf .Call()}. The
   former is inefficient and the latter requires knowledge of R
   internals.

   \item {\bf Rcpp} makes it easy.

   \end{itemize}

\item {\bf bigmemory}

   \begin{itemize}

   \item R currently not completely 64-bit. 

   \item Can have 52-bit integers, but only 32-bit matrix row/col
   dimensions.

   \item The {\bf bigmemory} package allows storing R matrices in ``C
   land,'' circumventing R storage limits.

   \item Storage is in {\bf shmem}, thus allowing for multicore use {\bf
   Rdsm}).

   \end{itemize}

\end{itemize}

\end{frame}

% \begin{frame}
% \frametitle{Rcpp Example}
% \pause
% 
% \end{frame}

\begin{frame}
\frametitle{Software Alchemy}
\pause

\begin{itemize}

\item For ``statistical'' problems, in ``iid'' form.
\pause
Image, text classification work.
\pause

\item Simple idea:

   \begin{itemize}

   \item Break data into ``independent'' chunks.

   \item Apply the procedure, e.g.\ logistic regression, to each chunk.

   \item Use combining op, e.g.\ averaging, for final answer.

   \item Provably correct and efficient.
   \pause

   \end{itemize}

\item A variant:  Apply procedure to chunks but take combining op to be
concatenation them rather than averaging.

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Serial Benefits of Software Alchemy}\pause

\begin{itemize}

\item SA gives speedup even in serial case of task is $O(n^c)$ for $c >
1$

\item Use SA to address a common problem: Big data, small GPU memory.
\pause
Apply GPU to each chunk, serially, then run combining op.
\pause

\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Example:  NMF}\pause

\begin{itemize}

\item E.g. break rows or columsn into $m$ chunks.

\item Get approximation $WH$ for each one.

\item To predict new case:

   \begin{itemize}

   \item Get the $m$ predictions.

   \item Combine via voting.

   \end{itemize}

\end{itemize}

\end{frame}

\end{document}