[解決済み] 本当に "*apply "ファミリーはベクトル化されていないのか?
質問
Rの新規ユーザーに対して、よくこう言いますね。
apply
がベクトル化されていないことを、Patrick Burnsの
R インフェルノ
サークル 4
"と書いてある(引用)。
一般的な反射神経として、applyファミリーの関数を使うことがあります。 これは ベクトル化ではなく、ループハディングです。 . apply関数の定義にforループがあります。 を定義しています。lapply関数はループを隠しますが、実行時間は明示的なforループとほぼ同じになる傾向があります。 実行時間は明示的なforループとほぼ等しくなる傾向があります。
確かに、ざっと見たところ
apply
のソースコードを見てみると、ループがあることがわかります。
grep("for", capture.output(getAnywhere("apply")), value = TRUE)
## [1] " for (i in 1L:d2) {" " else for (i in 1L:d2) {"
ここまではいいのですが、次に
lapply
または
vapply
というように、全く異なる図式が見えてきます。
lapply
## function (X, FUN, ...)
## {
## FUN <- match.fun(FUN)
## if (!is.vector(X) || is.object(X))
## X <- as.list(X)
## .Internal(lapply(X, FUN))
## }
## <bytecode: 0x000000000284b618>
## <environment: namespace:base>
ということで、どうやらRはないようです。
for
ループは存在せず、C言語で書かれた内部関数を呼び出しているようです。
をざっと見てみると ウサギ 穴 は、ほとんど同じ画像を表示します
さらに
colMeans
関数を例にとると、これはベクトル化されていないことを非難されることはありませんでした。
colMeans
# function (x, na.rm = FALSE, dims = 1L)
# {
# if (is.data.frame(x))
# x <- as.matrix(x)
# if (!is.array(x) || length(dn <- dim(x)) < 2L)
# stop("'x' must be an array of at least two dimensions")
# if (dims < 1L || dims > length(dn) - 1L)
# stop("invalid 'dims'")
# n <- prod(dn[1L:dims])
# dn <- dn[-(1L:dims)]
# z <- if (is.complex(x))
# .Internal(colMeans(Re(x), n, prod(dn), na.rm)) + (0+1i) *
# .Internal(colMeans(Im(x), n, prod(dn), na.rm))
# else .Internal(colMeans(x, n, prod(dn), na.rm))
# if (length(dn) > 1L) {
# dim(z) <- dn
# dimnames(z) <- dimnames(x)[-(1L:dims)]
# }
# else names(z) <- dimnames(x)[[dims + 1]]
# z
# }
# <bytecode: 0x0000000008f89d20>
# <environment: namespace:base>
あれ?これもただ単に
.Internal(colMeans(...
を呼び出すだけです。
ウサギの穴
. とはどう違うのでしょうか?
.Internal(lapply(..
?
実は、簡単なベンチマークで
sapply
よりも悪くないことがわかります。
colMeans
よりも悪くなく
for
ループよりも優れています。
m <- as.data.frame(matrix(1:1e7, ncol = 1e5))
system.time(colMeans(m))
# user system elapsed
# 1.69 0.03 1.73
system.time(sapply(m, mean))
# user system elapsed
# 1.50 0.03 1.60
system.time(apply(m, 2, mean))
# user system elapsed
# 3.84 0.03 3.90
system.time(for(i in 1:ncol(m)) mean(m[, i]))
# user system elapsed
# 13.78 0.01 13.93
つまり、次のように言ってよいのでしょうか。
lapply
と
vapply
は実際にはベクトル化された
(と比較して
apply
である
for
を呼び出すループです。
lapply
を呼び出す)、そしてパトリック・バーンズが本当に言いたかったことは何だったのでしょうか?
どのように解決するのか?
まず、あなたの例では "data.frame" でテストを行っていますが、これは以下のように公正ではありません。
colMeans
,
apply
と
"[.data.frame"
はオーバーヘッドがあるため
system.time(as.matrix(m)) #called by `colMeans` and `apply`
# user system elapsed
# 1.03 0.00 1.05
system.time(for(i in 1:ncol(m)) m[, i]) #in the `for` loop
# user system elapsed
# 12.93 0.01 13.07
マトリックスでは、少し様子が違います。
mm = as.matrix(m)
system.time(colMeans(mm))
# user system elapsed
# 0.01 0.00 0.01
system.time(apply(mm, 2, mean))
# user system elapsed
# 1.48 0.03 1.53
system.time(for(i in 1:ncol(mm)) mean(mm[, i]))
# user system elapsed
# 1.22 0.00 1.21
質問の主な部分をおさらいしておくと、この2つの主な違いは
lapply
/
mapply
/などと、ストレートなRループがループ処理を行う場所です。Rolandが指摘するように、CとRのループはどちらも各反復でR関数を評価する必要があり、これは最もコストがかかるものです。本当に高速なC関数は、すべてをCで行うものなので、quot;vectorised"とはこのことなのでしょう。
リスト("list")の各要素の平均を求める例です。
(
EDIT 5月11日 '16
: 私は、quot;mean"を求める例は、R関数の反復評価とコンパイルされたコードの違いのための良いセットアップではないと信じています、 (1) Rの平均アルゴリズムの特殊性から、quot;数値"の上に単純な
sum(x) / length(x)
でテストする方がより理にかなっていること、そして、(2)
length(x) >> lengths(x)
. そこで、"mean" の例は最後に移動して、別のものに置き換えています)。
簡単な例として、それぞれの反対を見つけることを考えることができます。
length == 1
の各要素の反対を見つけることを考えます。
で
tmp.c
ファイルを作成します。
#include <R.h>
#define USE_RINTERNALS
#include <Rinternals.h>
#include <Rdefines.h>
/* call a C function inside another */
double oppC(double x) { return(ISNAN(x) ? NA_REAL : -x); }
SEXP sapply_oppC(SEXP x)
{
SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));
for(int i = 0; i < LENGTH(x); i++)
REAL(ans)[i] = oppC(REAL(VECTOR_ELT(x, i))[0]);
UNPROTECT(1);
return(ans);
}
/* call an R function inside a C function;
* will be used with 'f' as a closure and as a builtin */
SEXP sapply_oppR(SEXP x, SEXP f)
{
SEXP call = PROTECT(allocVector(LANGSXP, 2));
SETCAR(call, install(CHAR(STRING_ELT(f, 0))));
SEXP ans = PROTECT(allocVector(REALSXP, LENGTH(x)));
for(int i = 0; i < LENGTH(x); i++) {
SETCADR(call, VECTOR_ELT(x, i));
REAL(ans)[i] = REAL(eval(call, R_GlobalEnv))[0];
}
UNPROTECT(2);
return(ans);
}
そしてR側では
system("R CMD SHLIB /home/~/tmp.c")
dyn.load("/home/~/tmp.so")
をデータで表示します。
set.seed(007)
myls = rep_len(as.list(c(NA, runif(3))), 1e7)
#a closure wrapper of `-`
oppR = function(x) -x
for_oppR = compiler::cmpfun(function(x, f)
{
f = match.fun(f)
ans = numeric(length(x))
for(i in seq_along(x)) ans[[i]] = f(x[[i]])
return(ans)
})
ベンチマークを行う。
#call a C function iteratively
system.time({ sapplyC = .Call("sapply_oppC", myls) })
# user system elapsed
# 0.048 0.000 0.047
#evaluate an R closure iteratively
system.time({ sapplyRC = .Call("sapply_oppR", myls, "oppR") })
# user system elapsed
# 3.348 0.000 3.358
#evaluate an R builtin iteratively
system.time({ sapplyRCprim = .Call("sapply_oppR", myls, "-") })
# user system elapsed
# 0.652 0.000 0.653
#loop with a R closure
system.time({ forR = for_oppR(myls, "oppR") })
# user system elapsed
# 4.396 0.000 4.409
#loop with an R builtin
system.time({ forRprim = for_oppR(myls, "-") })
# user system elapsed
# 1.908 0.000 1.913
#for reference and testing
system.time({ sapplyR = unlist(lapply(myls, oppR)) })
# user system elapsed
# 7.080 0.068 7.170
system.time({ sapplyRprim = unlist(lapply(myls, `-`)) })
# user system elapsed
# 3.524 0.064 3.598
all.equal(sapplyR, sapplyRprim)
#[1] TRUE
all.equal(sapplyR, sapplyC)
#[1] TRUE
all.equal(sapplyR, sapplyRC)
#[1] TRUE
all.equal(sapplyR, sapplyRCprim)
#[1] TRUE
all.equal(sapplyR, forR)
#[1] TRUE
all.equal(sapplyR, forRprim)
#[1] TRUE
(平均値の求め方の原型を踏襲しています)。
#all computations in C
all_C = inline::cfunction(sig = c(R_ls = "list"), body = '
SEXP tmp, ans;
PROTECT(ans = allocVector(REALSXP, LENGTH(R_ls)));
double *ptmp, *pans = REAL(ans);
for(int i = 0; i < LENGTH(R_ls); i++) {
pans[i] = 0.0;
PROTECT(tmp = coerceVector(VECTOR_ELT(R_ls, i), REALSXP));
ptmp = REAL(tmp);
for(int j = 0; j < LENGTH(tmp); j++) pans[i] += ptmp[j];
pans[i] /= LENGTH(tmp);
UNPROTECT(1);
}
UNPROTECT(1);
return(ans);
')
#a very simple `lapply(x, mean)`
C_and_R = inline::cfunction(sig = c(R_ls = "list"), body = '
SEXP call, ans, ret;
PROTECT(call = allocList(2));
SET_TYPEOF(call, LANGSXP);
SETCAR(call, install("mean"));
PROTECT(ans = allocVector(VECSXP, LENGTH(R_ls)));
PROTECT(ret = allocVector(REALSXP, LENGTH(ans)));
for(int i = 0; i < LENGTH(R_ls); i++) {
SETCADR(call, VECTOR_ELT(R_ls, i));
SET_VECTOR_ELT(ans, i, eval(call, R_GlobalEnv));
}
double *pret = REAL(ret);
for(int i = 0; i < LENGTH(ans); i++) pret[i] = REAL(VECTOR_ELT(ans, i))[0];
UNPROTECT(3);
return(ret);
')
R_lapply = function(x) unlist(lapply(x, mean))
R_loop = function(x)
{
ans = numeric(length(x))
for(i in seq_along(x)) ans[i] = mean(x[[i]])
return(ans)
}
R_loopcmp = compiler::cmpfun(R_loop)
set.seed(007); myls = replicate(1e4, runif(1e3), simplify = FALSE)
all.equal(all_C(myls), C_and_R(myls))
#[1] TRUE
all.equal(all_C(myls), R_lapply(myls))
#[1] TRUE
all.equal(all_C(myls), R_loop(myls))
#[1] TRUE
all.equal(all_C(myls), R_loopcmp(myls))
#[1] TRUE
microbenchmark::microbenchmark(all_C(myls),
C_and_R(myls),
R_lapply(myls),
R_loop(myls),
R_loopcmp(myls),
times = 15)
#Unit: milliseconds
# expr min lq median uq max neval
# all_C(myls) 37.29183 38.19107 38.69359 39.58083 41.3861 15
# C_and_R(myls) 117.21457 123.22044 124.58148 130.85513 169.6822 15
# R_lapply(myls) 98.48009 103.80717 106.55519 109.54890 116.3150 15
# R_loop(myls) 122.40367 130.85061 132.61378 138.53664 178.5128 15
# R_loopcmp(myls) 105.63228 111.38340 112.16781 115.68909 128.1976 15
関連
-
[解決済み] for'ループでインデックスにアクセスする?
-
[解決済み] callとapplyの違いは何ですか?
-
[解決済み] JavaScriptでオブジェクトのキー/プロパティの数を効率的にカウントする方法
-
[解決済み] Bashでファイルの中身をループする
-
[解決済み] オブジェクトをメンバーとして、プレーンなJavaScriptオブジェクトをループさせる方法
-
[解決済み] グループ化関数(tapply、by、aggregate)と*applyファミリ
-
[解決済み] INNER JOINよりもCROSS APPLYを使用すべきなのはどのような場合ですか?
-
[解決済み] 情報を損なわずに因数を整数値に変換するには?
-
[解決済み] Pandasのmap、applymap、applyメソッドの違い
-
[解決済み] Rのapplyファミリーは構文上の砂糖以上のものなのか?
最新
-
nginxです。[emerg] 0.0.0.0:80 への bind() に失敗しました (98: アドレスは既に使用中です)
-
htmlページでギリシャ文字を使うには
-
ピュアhtml+cssでの要素読み込み効果
-
純粋なhtml + cssで五輪を実現するサンプルコード
-
ナビゲーションバー・ドロップダウンメニューのHTML+CSSサンプルコード
-
タイピング効果を実現するピュアhtml+css
-
htmlの選択ボックスのプレースホルダー作成に関する質問
-
html css3 伸縮しない 画像表示効果
-
トップナビゲーションバーメニュー作成用HTML+CSS
-
html+css 実装 サイバーパンク風ボタン
おすすめ
-
R LanguageError in hist.default() : 'x' は数値でなければなりません.
-
R plot.new() のエラー : 図形の余白が大きすぎる
-
SocketTimeoutExceptionです。読み込みがタイムアウトしました
-
Rによる系統的クラスタリング(階層)分析のグラフ形式の完全版
-
[解決済み] lm.fit(x,y,offset = offset, singular.ok,...) 0 非NAケースでboxcox式で計算するとエラーになる。
-
[解決済み] Rでデータフレームに行を追加する方法は?
-
[解決済み] R dataframeでNAの値をゼロに置き換えるには?
-
[解決済み] 関数のソースコードを見るにはどうしたらいいですか?
-
[解決済み] なぜ `[`] は `subset` よりも優れているのですか?
-
[解決済み] Rのapplyファミリーは構文上の砂糖以上のものなのか?