gputools against base R: measured timings

The machine

The CPU side is base R as most people run it: a multi-threaded BLAS on 8 cores, which is what %*%, lm, qr and solve use. dist, hclust and cor do not use the BLAS and run on one core.

Each time is the median of three runs after a warm-up run, in seconds. Data is in R’s memory at the start and the result back in R’s memory at the end, so the GPU times include copying both ways.

Results

gpuMatMult only computes in double precision, so it has no single column.

#> Warning in xinch(0.1): x log scale: xinch() is nonsense

Reading the results

Where the GPU wins. Computations that are large, regular and many times over the same data: Kendall’s correlation (base R’s is a direct O(n^2) loop per pair of variables), Pearson’s and Spearman’s (one matrix product of standardised columns on the GPU), distance matrices, Granger tests over many pairs of series, the combined gpuDistClust, which keeps the distance matrix on the GPU for the clustering, and solving and decomposing large matrices (gpuSolve, gpuPrcomp).

Where it does not.

• gpuMatMult, gpuLm: a modern multi-threaded BLAS is already fast at these, and the GPU pays for copying the data across and back. The GPU here is also a workstation card with slow double precision.
• gpuHclust alone: the clustering is a sequence of n small dependent steps, which leaves the GPU mostly idle; it pays off only together with the distance computation (gpuDistClust).
• Small data of any kind: see the next section.

Single or double precision. Consumer and workstation GPUs run double precision many times slower than single; data-centre GPUs do not. For statistics whose result is a correlation, a distance or a test statistic, single precision usually loses nothing that matters; for least squares on badly conditioned data, use double.

When the GPU pays off

Every GPU call has a fixed cost: copying the data to the card and back, starting the computation, and the first time a function is used in a session, compiling its GPU code. On small data that fixed cost is most of the time, and base R is faster. Each function has a size above which the GPU wins. The table gives, for each, the smallest size tested at which it beat base R in single precision, and at which it was at least 5 times faster; the chart shows the whole range.

Filled points are single precision, open points double.

Wide data: the strongest pairs

For thousands of variables or more, the whole correlation matrix is the problem: 50,000 variables have 1,250,000,000 pairs, and their matrix takes 20 GB. Usually only the strong pairs are wanted. gpuCorPairs() returns those, every pair at or above a threshold or each variable’s strongest partners, computing the matrix on the GPU a block at a time and never holding it whole.

The table times finding every pair with |r| >= 0.8 among variables of 100 observations (random data with a few hundred planted correlated groups): base R’s cor() and then which(), gpuCor() and then which(), and gpuCorPairs(); and, last, gpuCorPairs() finding each variable’s 10 strongest partners. Times are in seconds; “–” was not run, the matrix being too big to wait for.

Most of the time cor() and gpuCor() take here goes on the matrix itself, building it and searching it in R, which gpuCorPairs() never does. Use it whenever you want pairs, not the matrix.

Association scans

gpuAssoc() tests every column of one matrix against every column of another by linear regression with covariates: the scan of an eQTL study, SNPs against genes. It computes what Matrix eQTL’s linear model does (Shabalin 2012), so the comparison is with Matrix eQTL itself, version 2.4, running on the multi-threaded OpenBLAS of the machine above. Each scan has 500 samples and 5 covariates and reports p <= 1e-5; both programs reported the same associations wherever both ran.

Times are in seconds; “–” was not run (at four billion tests Matrix eQTL ran out of a 48 GB job’s memory). At 1,000,000,000 tests the GPU was 30 times as fast. A genome-wide scan, millions of SNPs, is a loop over chromosomes or blocks of SNPs, each call large enough to keep the GPU busy.

K-means

gpuKmeans() runs Lloyd’s algorithm, step for step as kmeans(algorithm = "Lloyd") does. Base R’s default, Hartigan-Wong, is also timed, since it is what most people run: it does more work per pass but needs fewer passes, and usually ends at a lower total within-cluster sum of squares from the same start. All runs start from the same centres (50 dimensions, at most 20 iterations); the last column is the GPU’s sum of squares relative to Hartigan-Wong’s. Times are in seconds; “–” was not run, base R taking many minutes.

The GPU pays off from about ten thousand points and a hundred clusters, and most with many clusters. Where Lloyd’s algorithm ends at a worse sum of squares, a few random starts (nstart) cost little on the GPU and usually close the gap.

Permutation tests

A permutation test repeats the whole computation for every shuffle of the samples, so a thousand permutations cost a thousand times the test: the case the GPU is best at. Both functions here also give Westfall and Young’s maxT adjustment for testing everything at once, which needs every permutation’s results together.

gpuCorPerm(), 1000 permutations of 50 samples, every column of x against every column of y, against the same test in R: a loop of cor() calls, and the faster loop of BLAS products of standardised columns. Times in seconds; “–” was not run.

gpuRowTtest(), Welch’s t test of every gene between two groups of samples with 1000 permutations, against the same test in vectorised R (each permutation’s t for every gene from rowMeans and rowSums). The t tests themselves, without permutations, run on the CPU and take the time in the last column.

Your own hardware

To measure your own machine, run from a checkout of the source repository with gputools installed:

Rscript bench/run-benchmarks.R my-results