Rcpp 套件：結合 R 與 C++ 加速程式執行

Rcpp 是一個可以讓 R 透過 C++ 語言提升執行效能的套件，其使用方式簡單、上手容易，且對於耗時的運算相當有幫助，是一個很實用的套件。

由於 R 屬於高階語言，相較於 C/C++ 這類的低階語言，R 程式的執行速度比較慢，對於較為耗時的運算（例如蒙地卡羅類型的模擬）來說，若只使用單純的 R 語言來處理，時常會遇到計算時間過長的問題。

傳統上運算速度的問題可以使用 R 本身所提供的 C 或 Fortran 語言介面，將程式中耗時的部分以 C 或 Fortran 改寫，然而 R 本身的低階語言介面相當複雜，對於一般人而言不容易使用，通常只有 R 的套件開發者有能力使用這樣的功能。

Rcpp 是後來發展出來的一個套件，它在 R 中提供了一個非常簡單又好用的 C++ 語言介面，讓一般的使用者可以很輕鬆的結合 R 與 C++ 程式，降低程式開發的技術門檻，大幅提高 R 程式的執行速度，近年來 Rcpp 已經成為 R 高效能運算領域非常重要的套件之一。

適合使用 C++ 的問題

C++ 語言雖然執行效率比 R 語言更好，不過並不是每一種問題都適合使用 C++ 來處理，以下是幾種非常適合使用 C++ 來處理的典型問題類型：

無法轉換成向量化的迴圈運算，例如每次迭代都會跟前一次迭代的運算結果有關的迴圈。
遞迴呼叫的函數或是需要大量重複執行的函數，因為在 C++ 中的函數呼叫會比 R 中的函數呼叫更有效率，所以這種狀況也適合以 C++ 改寫。
需要使用 R 所沒有支援的複雜資料結構時，也可以使用 C++ 與其標準函式庫（STL）來實做。

安裝

由於 Rcpp 套件在使用時會需要編譯 C++ 的程式碼，所以除了安裝 R 的 Rcpp 套件之外，還要再安裝系統上的 C++ 編譯器才能正常使用。

安裝 Rtools 工具集

由於在 Windows 系統上通常並不會有安裝 C++ 編譯器，若要使用 Rcpp 套件，必須先安裝 Rtools 這個官方的工具集，而若是在 Linux 系統上的話，請安裝一般的 gcc/g++ 編譯器。

Rtools 工具集裡面包含了編譯 R 套件所需要的各種工具，而我們需要的就是其中的 gcc/g++ 編譯器。

Step 1

請從官方的 CRAN 網站上下載 Rtools 安裝檔。

r-install-rtools-1

Step 2

下載時請選擇適合自己 R 版本的 Rtools 安裝檔。

r-install-rtools-2

Step 3

下載下來之後，直接執行，進行安裝。首先選擇安裝用語言，選擇 English 即可，按下「OK」。

r-install-rtools-3

Step 4

這是 Rtools 工具集的說明資訊，其中最主要的就是 Cygwin 的仿 Linux 環境，請點選「Next」繼續。

r-install-rtools-4

Step 5

選擇安裝路徑，點選「Next」繼續。

r-install-rtools-5

Step 6

選擇安裝組件，使用預設值即可，點選「Next」繼續。

r-install-rtools-6

Step 7

其他設定，使用預設值即可，點選「Next」繼續。

r-install-rtools-7

Step 8

檢查安裝設定，確認無誤後，點選「Install」進行安裝。

r-install-rtools-8

Step 9

等待安裝過程。

r-install-rtools-9

Step 10

安裝完成，點選「Finish」結束。

r-install-rtools-10

安裝 Rcpp 套件

R 官方的套件庫就有收錄 Rcpp 套件，直接以 install.packages 安裝即可：

install.packages("Rcpp")

安裝完成後，即可載入使用：

library(Rcpp)

開始使用 Rcpp

Rcpp 的 cppFunction 是一個可以讓您將 C++ 程式碼直接內崁至 R 程式碼中的包裝函數，其用法非常簡單，在詳細介紹之前，我們先來看一個典型的範例，先讓大家對 Rcpp 的使用有一個初步的概念。

首先將 R 程式中需要大量運算的部分，改以 C++ 語言來撰寫：

int cpp_add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}

這個 C++ 函數是傳入 x、y 與 z 三個整數，然後計算它們的總和（sum），最後將計算結果傳回。

撰寫好 C++ 的計算函數之後，直接將 C++ 的整個程式碼內容當成參數，傳遞給 R 的 cppFunction 函數：

cppFunction('int cpp_add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')

cppFunction 會對這段 C++ 程式碼自動進行編譯，並且載入至目前的 R 環境中，若 C++ 的程式碼沒有錯誤的話，這時候在 R 中就會建立好一個剛剛 C++ 中所定義的函數：

cpp_add

function (x, y, z)
.Primitive(".Call")(, x, y, z)

而這個時候我們就可以在 R 中直接使用這個以 C++ 撰寫的 cpp_add 函數了：

cpp_add(1, 2, 3)

[1] 6

以上就是一個簡單的 Rcpp 使用範例，從這個範例您可以看得出來 Rcpp 的使用方式相當簡潔，不需要任何手動編譯的過程。

由於使用 Rcpp 的同時也必須撰寫一些 C++ 程式，若事先對於 C++ 語言有基本的認識會有很大的幫助。

以下我們將詳細介紹各種參數類型的用法，包含各類型的純量、向量、矩陣的輸入與輸出方式。

無輸入、純量輸出

一個沒有輸入值、只有單一傳回值的 R 函數是這樣寫的：

one <- function() 1L

而若將這個 R 函數以 C++ 改寫，則會像這樣：

int one() {
  return 1;
}

這裡的 int（integer）是指定此函數的傳回值型態為整數，而 return 則是指定實際要傳回的數值。

在 R 函數中可以省略 return 這一行，以最後一個運算結果作為函數傳回值，但是在 C++ 中則不可省略，一定要使用 return 明確指定函數的傳回值。另外 C++ 中的每一行運算式結尾都一定要加上分號（;），不可省略。

在 R 中我們可以將這個 C++ 函數透過 cppFunction 編譯並載入使用：

cppFunction('int one() {
  return 1;
}')
one()

[1] 1

在使用 cppFunction 編譯 C++ 函數時，它會自動在 R 中建立一個相同名稱的 R 函數，以 .Call 呼叫這個編譯好的 C++ 函數，使用者不需要自行建立這個 R 函數。

C++ 函數在定義時需要明確指定其傳回值的型態，若傳回值為簡單的純量（scalar），可以使用 C++ 內建的資料型別：

C++ 語法	說明
`double`	浮點數
`int`	整數
`std::string`	字串
`bool`	布林值

例如一個傳回字串的 C++ 函數就可以這樣寫：

std::string my_name() {
  return "G. T. Wang";
}

使用 cppFunction 編譯並載入使用：

cppFunction('std::string my_name() {
  return "G. T. Wang";
}')
my_name()

[1] "G. T. Wang"

純量輸入、純量輸出

下面這個 R 函數是一個以勾股定理計算直角三角形斜邊長度的函數，輸入三角形的勾長與股長，傳回弦長（斜邊的長度）：

hypotenuse <- function(a, b) {
  sqrt(a * a + b * b)
}

若以 C++ 改寫，則為：

double hypotenuse(double a, double b) {
  return sqrt(a * a + b * b);
}

C++ 函數的輸入參數與傳回值一樣都需要明確指定變數的類型，由於在計算斜邊長度時，會需要浮點運算，所以這裡我們將輸入與輸出的變數都宣告為浮點數（double），而輸入參數可使用的變數類型跟傳回值可使用的變數類型相同（double、int 等）。

這裡我們用到一個 C++ 的平方根函數 sqrt，其使用方式跟 R 中的 sqrt 相同，可計算傳入數值的平方根。

使用 cppFunction 編譯並載入使用：

cppFunction('double hypotenuse(double a, double b) {
  return sqrt(a * a + b * b);
}')
hypotenuse(3, 4)

[1] 5

向量輸入、純量輸出

C++ 語言的一個很大的優勢就是它的迴圈計算速度比 R 快很多，一個典型的用法是將 R 中的向量傳入 C++ 中進行運算，最後再傳回運算的結果。下面這個例子是一個計算總和的 R 函數：

r_sum <- function(x) {
  total <- 0
  for (i in seq_along(x)) {
    total <- total + x[i]
  }
  total
}

若要將這個 R 函數以 C++ 改寫，則會變成：

double cpp_sum(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}

由於這裡傳入的參數 x 會是一個 R 的向量，在 C++ 中純量與向量是不同的資料型別，Rcpp 針對 R 的幾種向量特別定義了一些 C++ 的類別供使用者使用，這裡我們將 x 的型別宣告為浮點數向量 NumericVector，並且使用 NumericVector 的 size 方法函數取得向量的長度，接著宣告一個用來儲存總和值的 total 變數（由於我們要計算浮點數向量的總和，因此將 total 變數宣告為浮點數 double），然後進行後續的 for 迴圈運算，最後將計算完成的總和值傳回。

以下是各種 R 向量所對應的 C++ 類別：

C++ 類別	說明
`NumericVector`	浮點數向量
`IntegerVector`	整數向量
`CharacterVector`	字元向量
`LogicalVector`	邏輯向量

使用 cppFunction 編譯並載入使用：

cppFunction('double cpp_sum(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}')

我們可以使用 R 的 microbenchmark 套件來測試 R 與 C++ 之間的迴圈運算效能差異：

library(microbenchmark)
x <- runif(1e3)
x.benchmark <- microbenchmark(
  sum(x),
  cpp_sum(x),
  r_sum(x)
)
x.benchmark

Unit: microseconds
       expr     min       lq      mean   median       uq      max neval
     sum(x)   1.131   1.3255   1.61292   1.5745   1.7285    3.708   100
 cpp_sum(x)   2.832   3.2555   4.20410   3.9425   4.7900   14.222   100
   r_sum(x) 344.723 389.2560 444.48158 414.5810 442.0710 1255.476   100

這裡我們測試了一般的 R 函數（r_sum）、C++ 函數（cpp_sum）以及 R 內建的向量化函數（sum），內建的向量化函數執行效率是最好的，而我們自己定義的 C++ 函數稍微遜色一些，而普通 R 函數的執行速度則是比我們定義的 C++ 函數慢了百倍以上。

畫出測試結果的小提琴圖（violin plot）：

library(ggplot2)
autoplot(x.benchmark)

測試結果小提琴圖

向量輸入、向量輸出

輸入與輸出都是向量的函數也很常見，這裡我們以一個氣泡排序演算法的排序函數做示範，以下是一個 R 版本的排序函數：

r_sort <- function(x) {
  x.length = length(x)
  for (i in seq(1, x.length-1)) {
    for (j in seq(1, x.length-i)) {
      if (x[j] > x[j+1]) {
        tmp = x[j]
        x[j] = x[j+1]
        x[j+1] = tmp
      }
    }
  }
  x
}

同樣的演算法以 C++ 改寫則變成：

NumericVector cpp_sort(NumericVector x) {
  int x_length = x.size();
  for(int i = 0; i < x_length-1; ++i) {
    for(int j = 0; j < x_length-1-i; ++j) {
      if (x[j] > x[j+1]) {
        double tmp = x[j];
        x[j] = x[j+1];
        x[j+1] = tmp;
      }
    }
  }
  return x;
}

使用 cppFunction 編譯並載入使用：

cppFunction('NumericVector cpp_sort(NumericVector x) {
  int x_length = x.size();
  for(int i = 0; i < x_length-1; ++i) {
    for(int j = 0; j < x_length-1-i; ++j) {
      if (x[j] > x[j+1]) {
        double tmp = x[j];
        x[j] = x[j+1];
        x[j+1] = tmp;
      }
    }
  }
  return x;
}')

在實作這種比較複雜一點的演算法時，通常在程式寫完之後，最好先做一下驗證，確保程式沒有寫錯。這裡我們將自己寫的兩個排序函數跟 R 內建的 sort 函數比較，看看結果是否完全相同：

x <- runif(100)
x.sort <- sort(x)
x.r_sort <- r_sort(x)
x.cpp_sort <- cpp_sort(x)
sum((x.sort - x.r_sort)^2)

[1] 0

sum((x.sort - x.cpp_sort)^2)

[1] 0

確認程式無誤後，使用 microbenchmark 測試執行效能：

library(microbenchmark)
x <- runif(100)
x.benchmark <- microbenchmark(
  sort(x),
  cpp_sort(x),
  r_sort(x)
)
x.benchmark

Unit: microseconds
        expr      min        lq       mean    median        uq      max neval
     sort(x)   11.933   13.1045   21.55053   19.1545   31.8010   44.421   100
 cpp_sort(x)    6.723    7.3220    9.91120    9.7080   12.6855   20.863   100
   r_sort(x) 3626.665 3713.8905 3972.78264 3816.2110 3920.3770 5933.124   100

畫出測試結果的小提琴圖（violin plot）：

library(ggplot2)
autoplot(x.benchmark)

測試結果小提琴圖

在這個例子中，自己使用 Rcpp 寫的 C++ 排序函數，其執行效能還比內建的 sort 函數還要好，而一般的 R 排序函數執行速度一樣是最慢的。若將 x 向量的長度增加一些，結果又會有些不同：

library(microbenchmark)
x <- runif(1000)
x.benchmark <- microbenchmark(
  sort(x),
  cpp_sort(x),
  r_sort(x)
)
x.benchmark

Unit: microseconds
        expr        min          lq        mean      median          uq        max neval
     sort(x)     16.569     18.0955     31.3704     40.2485     43.1295     54.350   100
 cpp_sort(x)    380.192    381.5215    391.0866    388.8055    390.3740    648.092   100
   r_sort(x) 302941.105 304085.3015 314353.2839 304802.2955 306465.2875 751781.425   100

library(ggplot2)
autoplot(x.benchmark)

測試結果小提琴圖

對於長度比較長的 x 向量，內建的 sort 排序函數就會比較快一些，而一般的 R 排序函數執行速度還是一如往常的慢。

氣泡排序是比較簡單的演算法，執行所需時間為 O((n^2))，而內建的 sort 函數所使用的排序演算法絕對比 O((n^2)) 還要好，所以在向量長度較長時，差異會比較明顯。

矩陣輸入、向量輸出

Rcpp 同時也對 R 的幾種矩陣定義了一些 C++ 的類別：

C++ 類別	說明
`NumericMatrix`	浮點數矩陣
`IntegerMatrix`	整數矩陣
`CharacterMatrix`	字元矩陣
`LogicalMatrix`	邏輯矩陣

其使用方式跟向量類似，以下是將 R 的 rowSums 改寫成 C++ 的程式碼：

NumericVector cpp_rowSums(NumericMatrix x) {
  int nrow = x.nrow(), ncol = x.ncol();
  NumericVector result(nrow);
  for (int i = 0; i < nrow; i++) {
    double total = 0;
    for (int j = 0; j < ncol; j++) {
      total += x(i, j);
    }
    result[i] = total;
  }
  return result;
}

NumericMatrix 的 nrow 與 ncol 方法可以用來取得該矩陣的列（row）數與行（column）數，而在提取矩陣內部的元素時，是使用小括號 () 運算子。

使用 cppFunction 編譯並載入使用：

cppFunction('NumericVector rowSumsC(NumericMatrix x) {
  int nrow = x.nrow(), ncol = x.ncol();
  NumericVector result(nrow);
  for (int i = 0; i < nrow; i++) {
    double total = 0;
    for (int j = 0; j < ncol; j++) {
      total += x(i, j);
    }
    result[i] = total;
  }
  return result;
}')

與 R 內建的 rowSums 函數做比較：

set.seed(1)
x <- matrix(sample(100), 10)
rowSums(x)

[1] 524 494 436 529 423 505 553 451 577 558

cpp_rowSums(x)

[1] 524 494 436 529 423 505 553 451 577 558

`sourceCpp` 函數

前面我們介紹過以 cppFunction 函數將 C++ 程式碼內崁在 R 程式碼中，這樣的方式對於簡短的小程式而言非常方便，但是如果需要開發比較複雜的 C++ 程式時，把大量的 C++ 程式碼貼在 cppFunction 函數中，會讓程式碼難以閱讀，在開發上也會很不方便。

遇到比較複雜的程式架構時，我們可以將 C++ 程式碼獨立出來，另外儲存成一個 .cpp 檔案，然後使用 sourceCpp 函數來編譯並載入 .cpp 檔案中的 C++ 程式碼，這樣的作用跟 cppFunction 函數相同，但是可以讓大型的程式更好管理。

當我們將 C++ 的程式碼獨立出來時，必須引入 Rcpp.h 並且設定 Rcpp 這個 namespace，也就是在檔案的開頭加上這兩行：

#include <Rcpp.h>
using namespace Rcpp;

對於那些需要在 R 中被呼叫的 C++ 函數之前，加上這一行註解：

// [[Rcpp::export]]

請注意這一行註解中包含的空白不可以省略。

我們將前面計算總和的 C++ 範例程式改成獨立的一個 cpp_sum.cpp 檔，內容如下：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double cpp_sum(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}

接著以 sourceCpp 編譯並載入使用：

sourceCpp("cpp_sum.cpp")
cpp_sum(1:10)

[1] 55

我們也可以將 R 的程式碼內崁在 C++ 的註解之中，其語法為：

/*** R
# 這是 R 的程式碼
*/

放在這裡的 R 程式碼會以 source(echo = TRUE) 的方式來執行，所以我們不需要特別以 print 等函數來輸出，這個功能對於開發階段的測試與除錯會有一些幫助。以下是計算總和的範例：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double cpp_sum(NumericVector x) {
  int n = x.size();
  double total = ;
  for(int i = ; i < n; ++i) {
    total += x[i];
  }
  return total;
}

/*** R
library(microbenchmark)
x <- runif(1e3)
microbenchmark(
  sum(x),
  cpp_sum(x)
)
*/

將這段 C++ 程式碼儲存在 cpp_sum_embedded_r.cpp 檔案中，然後同樣以 sourceCpp 編譯並載入：

sourceCpp("cpp_sum_embedded_r.cpp")

> library(microbenchmark)
> x <- runif(1e3)
> microbenchmark(
+   sum(x),
+   cpp_sum(x)
+ )
Unit: microseconds
       expr   min    lq    mean median     uq    max neval
     sum(x) 1.136 1.161 1.53268 1.2245 1.2485 29.916   100
 cpp_sum(x) 2.689 2.738 3.12136 2.7840 2.8880 23.875   100

當 C++ 的函數被載入之後，內崁的 R 程式碼也會一併被執行。

範例程式碼

這裡對於一些常見的程式結構提供了 R 與 C++ 的對照範例程式碼，對於 C++ 不熟新的使用者可以參考使用。

`if` 判斷式

這是使用 if 判斷式計算絕對值的 R 函數：

r_abs <- function(x) {
  if (x > ) {
    x
  } else {
    -x
  }
}

C++ 版本則為：

double cpp_abs(double x) {
  if (x > 0) {
    return x;
  } else {
    return -x;
  }
}

`switch` 判斷式

這是一段 R 語言的 switch 判斷式範例程式碼：

r_int2str <- function(x) {
  switch(as.character(x),
    "1" = "one",
    "2" = "two",
    "oops"
  )
}

C++ 版本則為：

std::string cpp_int2str(int x) {
  switch(x) {
    case 1:
      return "one";
    case 2:
      return "two";
    default:
      return "oops";
  }
}

`for` 迴圈

這是使用 for 迴圈計算平均值的 R 函數：

r_mean <- function(x) {
  n = length(x)
  total = 0
  for(i in 1:n) {
    total <- total + x[i]
  }
  total / n
}

C++ 版本則為：

double cpp_mean(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total / n;
}

`while` 迴圈

這是使用 while 迴圈計算階乘的 R 函數：

r_factorial <- function(x) {
  result = x
  while (x > 2) {
    x <- x - 1
    result <- result * x
  }
  result
}

C++ 版本則為：

int cpp_factorial(int x) {
  int result = x;
  while (x > 2) {
    x -= 1;
    result *= x;
  }
  return result;
}

`next` 與 `break`

這是使用 next 與 break 列出奇數的 R 函數：

r_odd <- function(x) {
  i <- 0
  result = integer()
  while ( TRUE ) {
    i <- i + 1
    if (i > x) break
    if (i %% 2 == 0) next
    result = c(result, i)
  }
  result
}

C++ 版本則為：

IntegerVector cpp_odd(int x) {
  int i = 0;
  IntegerVector result;
  while ( true ) {
    i += 1;
    if (i > x) break;
    if (i % 2 == 0) continue;
    result.push_back(i);
  }
  return result;
}

IntegerVector 的 push_back 跟 std::vector 的 push_back 用法相同，可將新元素從後方加入向量中。

🚧 文件撰寫中…

物件屬性

列表（List）與 data frames

函數

其他類型物件

缺失值

純量

字串

布林值

向量

Rcpp sugar

參考資料

Advanced R

https://cran.rstudio.com/web/packages/Rcpp/vignettes/Rcpp-introduction.pdf https://support.rstudio.com/hc/en-us/articles/200486088-Using-Rcpp-with-RStudio https://cran.rstudio.com/web/packages/Rcpp/vignettes/Rcpp-quickref.pdf https://cran.rstudio.com/web/packages/Rcpp/vignettes/Rcpp-attributes.pdf http://gallery.rcpp.org/ http://web.ydu.edu.tw/~alan9956/R/Rcpp-Rtools-tutorial.pdf

適合使用 C++ 的問題#

安裝#

安裝 Rtools 工具集#

安裝 Rcpp 套件#

開始使用 Rcpp#

無輸入、純量輸出#

純量輸入、純量輸出#

向量輸入、純量輸出#

向量輸入、向量輸出#

矩陣輸入、向量輸出#

sourceCpp 函數#

範例程式碼#

if 判斷式#

switch 判斷式#

for 迴圈#

while 迴圈#

next 與 break#

物件屬性#

列表（List）與 data frames#

函數#

其他類型物件#

缺失值#

純量#

字串#

布林值#

向量#

Rcpp sugar#

參考資料#

適合使用 C++ 的問題

安裝

安裝 Rtools 工具集

安裝 Rcpp 套件

開始使用 Rcpp

無輸入、純量輸出

純量輸入、純量輸出

向量輸入、純量輸出

向量輸入、向量輸出

矩陣輸入、向量輸出

`sourceCpp` 函數

範例程式碼

`if` 判斷式

`switch` 判斷式

`for` 迴圈

`while` 迴圈

`next` 與 `break`

物件屬性

列表（List）與 data frames

函數

其他類型物件

缺失值

純量

字串

布林值

向量

Rcpp sugar

參考資料