从R调用Julia性能极差的原因排查求助

测试脚本 #

# Loading packages
library(Rcpp)
library(RcppArmadillo)
library(JuliaConnectoR)
library(microbenchmark)

# Generating data
set.seed(9392927)
N <- 1e3L
x1 <- matrix(rpois(N ^ 2L, lambda = 10L), N)
x2 <- matrix(rpois(N ^ 2L, lambda = 10L), N)

# Defining functions
mySolveR <- function(m) {
  return(solve(m))
}
myProdR <- function(m1, m2) {
  return(m1 %*% m2)
}
cppFunction("
  arma::mat mySolveCpp(arma::mat m) {
    return(inv(m));
  }",
  depends = "RcppArmadillo"
)
cppFunction("
  arma::mat myProdCpp(arma::mat m1, arma::mat m2) {
    return(m1 * m2);
  }",
  depends = "RcppArmadillo"
)
mySolveJulia <- juliaEval("
  function mySolveJulia(m)
    return(inv(m))
  end"
)
#> Starting Julia ...
myProdJulia <- juliaEval("
  function myProdJulia(m1, m2)
    return(m1 * m2)
  end"
)

# Benchmarking
reps <- 10L
print(
  microbenchmark(
      mySolveR(x1), mySolveCpp(x1), mySolveJulia(x1),
      times = reps
    )
)
#> Unit: milliseconds
#>              expr        min         lq      mean     median        uq
#>      mySolveR(x1)   41.53371   44.03872  186.4527   59.43063  465.0983
#>    mySolveCpp(x1)   49.42867   72.63049  240.0839   81.16933  440.1176
#>  mySolveJulia(x1) 3856.43147 3897.36013 5112.4448 5376.94133 6026.1369
#>        max neval cld
#>   528.8838    10  a 
#>   576.2336    10  a 
#>  6220.1093    10   b
print(
  microbenchmark(
    myProdR(x1, x2), myProdCpp(x1, x2), myProdJulia(x1, x2),
    times = reps
  )
)
#> Unit: milliseconds
#>                 expr        min         lq       mean     median         uq
#>      myProdR(x1, x2)   20.10827   23.50031   27.69398   25.18961   32.51614
#>    myProdCpp(x1, x2)   25.59671   29.59619   38.32624   31.83016   49.33401
#>  myProdJulia(x1, x2) 5698.10212 6135.62034 6353.71632 6412.39110 6487.12516
#>         max neval cld
#>    39.89406    10  a 
#>    57.82850    10  a 
#>  7243.81362    10   b

问题描述 #

原本预期Julia在矩阵求逆、矩阵乘法这类基础运算上的性能应与R、C处于同一数量级，但测试结果显示Julia的耗时是R/C的几十倍。已尝试startJuliaServer()/stopJulia()、将函数单独放在源文件等方法，问题仍未解决。

核心原因与解决方案 #

1. 未进行JIT预热编译 #

Julia的JIT编译器会在函数第一次调用时完成编译，这个编译时间会被计入第一次调用的耗时。基准测试直接开始计时，前几次调用的耗时包含了编译开销，拉高了整体平均值。

解决方法：在基准测试前先手动调用一次Julia函数，完成预热：

# 预热JIT编译
mySolveJulia(x1)
myProdJulia(x1, x2)

# 再执行基准测试
print(microbenchmark(...))

2. 重复的数据序列化/传输开销 #

每次调用mySolveJulia(x1)或myProdJulia(x1,x2)时，R都需要将1000×1000的大矩阵序列化后传输到Julia进程，计算完成后再将结果序列化传回R。这部分开销远大于计算本身，是导致Julia耗时过高的主要原因。

解决方法：将数据一次性传输到Julia环境，后续调用直接操作Julia中的数据：

# 将R中的矩阵赋值到Julia环境
juliaAssign("x1_jl", x1)
juliaAssign("x2_jl", x2)

# 定义操作Julia本地数据的函数
mySolveJulia_local <- juliaEval("
  function mySolveJulia_local()
    return(inv(x1_jl))
  end
")
myProdJulia_local <- juliaEval("
  function myProdJulia_local()
    return(x1_jl * x2_jl)
  end
")

# 预热后测试
mySolveJulia_local()
myProdJulia_local()

print(
  microbenchmark(
      mySolveR(x1), mySolveCpp(x1), mySolveJulia_local(),
      times = reps
    )
)
print(
  microbenchmark(
    myProdR(x1, x2), myProdCpp(x1, x2), myProdJulia_local(),
    times = reps
  )
)

3. 函数类型稳定性不足 #

R传入的矩阵类型在Julia中是Matrix{Int64}，而inv()和矩阵乘法需要浮点数矩阵，Julia会自动进行类型转换，但如果函数没有明确的类型声明，可能会影响JIT编译的优化效果。

解决方法：给Julia函数添加明确的类型注解，提前完成类型转换：

mySolveJulia <- juliaEval("
  function mySolveJulia(m::Matrix{Int64})
    return(inv(convert(Matrix{Float64}, m)))
  end
")
myProdJulia <- juliaEval("
  function myProdJulia(m1::Matrix{Int64}, m2::Matrix{Int64})
    return(convert(Matrix{Float64}, m1) * convert(Matrix{Float64}, m2))
  end
")

通过以上调整，Julia的计算性能应该能和R/C++处于同一数量级，体现出其本身的运算效率。

内容的提问来源于stack exchange，提问作者Waldir Leoncio