この一つ前の記事で。
二重の for は効率が悪い。outer を使うのがよい。
system.time({
x <- seq(-0.5, 3, length = 1000)
y <- seq(-2.5, 2.5, length = 1000)
z <- matrix(0, 1000, 1000)
for (i in 1:1000) {
for (j in 1:1000) {
z[i, j] <- (x[i] ^ 2 + y[j] ^ 2 - 2 * x[i]) ^ 3 - 4 * x[i] ^ 2 * y[j] ^ 2
}
}
contour(x, y, z, drawlabels = FALSE, levels = c(0, 1, 2, 5, 10, 20, 30, 50))
})
# ユーザ システム 経過
# 11.962 0.119 12.111
10 倍速だ。
system.time({
x <- seq(-0.5, 3, length = 1000)
y <- seq(-2.5, 2.5, length = 1000)
z <- outer(x, y, function(xi, yj) (xi ^ 2 + yj ^ 2 - 2 * xi) ^ 3 - 4 * xi ^ 2 * yj ^ 2)
contour(x, y, z, drawlabels = FALSE, levels = c(0, 1, 2, 5, 10, 20, 30, 50))
})
# ユーザ システム 経過
# 0.923 0.037 0.956
下記のような cxxfunction を使うのと遜色ない実行速度
library(inline)
src <- '
Rcpp::NumericVector x(X);
Rcpp::NumericVector y(Y);
int nx = x.size();
int ny = y.size();
Rcpp::NumericMatrix z(nx, ny);
for (int i = 0; i < nx; i++) {
for (int j = 0; j < ny; j++) {
z(i, j) = pow(x(i) * x(i) + y(j) * y(j) - 2 * x(i), 3) - 4 * x(i) * x(i) * y(j) * y(j);
}
}
return wrap(z);
'
calcZ <- cxxfunction(signature(X="numeric", Y="numeric"), src, plugin="Rcpp")
x <- seq(-0.5, 3, length = 1000)
y <- seq(-2.5, 2.5, length = 1000)
z <- calcZ(x, y)
contour(x, y, z, drawlabels = FALSE, levels = c(0, 1, 2, 5, 10, 20, 30, 50))
もっとも，これを以下のようにするとちょっとだけ速くなる。ループ内で不変な計算はループの外でやることと，何回も同じ計算をするときには一回計算して変数に代入しておくこと。
src <- '
Rcpp::NumericVector x(X);
Rcpp::NumericVector y(Y);
int nx = x.size();
int ny = y.size();
Rcpp::NumericMatrix z(nx, ny);
for (int i = 0; i < nx; i++) {
double xi2 = x(i) * x(i);
for (int j = 0; j < ny; j++) {
double yj2 = y(j) * y(j);
z(i, j) = pow(xi2 + yj2 - 2 * x(i), 3) - 4 * xi2 * yj2;
}
}
return wrap(z);
'
