Rもだいたい線形代数の範囲で片付けている限りは、そこそこ速いことがわかります。Fortranとの速度差が出ないというか、R vs R/FortranだとRの方が微妙に速く*1。
遅くなるのになぜFortranでRの対数尤度関数とその一階微分を用意したかと言うと、Fortranでハミルトニアン・モンテカルロ法を試みるために必要なコードで、準備として整備する必要があるからです。中間段階の副産物として、大してループが無い行列の計算でまとまるものは、Rで書いておけば済むことがわかったので記します*2。
正規分布の対数尤度関数とその解析的な一階微分
BLASを使うか否かを迷ったのですが、非BLASの方が速かったので非BLAS版を使いました。
module norm_d implicit none contains double precision function llf_norm(nr, nc, X, y, p) implicit none integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p integer :: i, j double precision, dimension(nr) :: s !$omp parallel !$omp do do i = 1, nr s(i) = dlog(dnorm(y(i) - sum(X(i, :)*p(2: 1 + nc)),p(1))) end do !$omp end do !$omp end parallel llf_norm = sum(s) end function double precision function dnorm(x, sigma) implicit none double precision, intent(in) :: x, sigma double precision, parameter :: pi = 3.141592653589793115998 dnorm = dexp(-x**2/2/sigma**2)/dsqrt(2*pi*sigma**2) end function subroutine llfg_norm(nr, nc, X, y, p, g) implicit none integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p double precision, dimension(nc + 1), intent(out) :: g double precision, dimension(nc + 1) :: q double precision :: theta, dtheta double precision, dimension(nc) :: beta, dbeta double precision, dimension(nr) :: res integer :: i beta(:) = p(2:nc + 1) theta = p(1) ! BLAS ! res(:) = y(:) ! call dgemm('n', 'n', nr, 1, nc, -1d0, X, nr, beta, nc, 1d0, res, nr) ! No BLAS !$omp parallel !$omp do do i = 1, nr res(i) = y(i) - sum(X(i, :)*beta) end do !$omp end do !$omp end parallel dtheta = -1*dble(nr)/theta + sum(res(:)**2)/theta**3 do i = 1, nc dbeta(i) = sum(X(:, i)*2*res)/theta**2/2 end do g(1) = dtheta g(2:nc + 1) = dbeta end subroutine end module
ファイル名を llf_norm.f90 として保存します。
数値微分のコード
一階微分のコードを書いて検算をするのは手間なので、こちらの数値微分のコードを使う予定ですが、上の解析的な微分のコードを使う方が1.5倍から2倍ほど速かったです。
subroutine diff_llf(llf, nr, nc, X, y, p, g) implicit none interface double precision function llf(nr, nc, X, y, p) implicit none integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p end function end interface integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p double precision, dimension(nc + 1), intent(out) :: g double precision, dimension(nc + 1) :: q double precision :: h = 1d-6, v_h, v_l integer :: i q(:) = p(:) ! 入力データはいじらない do i = 1, nc + 1 q(i) = p(i) + h v_h = llf(nr, nc, X, y, q) q(i) = p(i) - h v_l = llf(nr, nc, X, y, q) g(i) = (v_h - v_l)/(2*h) q(i) = p(i) ! 戻す end do end subroutine
ファイル名を diff_llf.f90 として保存します。
Rから呼び出すFortranのプロシージャ
Fortranの関数やプロシージャを何でも呼べるわけではないので、Rから呼び出すFortranのプロシージャを用意します。
subroutine llf(r, nr, nc, X, y, p) use norm_d, only:llf_norm implicit none double precision, intent(out) :: r integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p r = llf_norm(nr, nc, X, y, p) end subroutine subroutine llfg(r, nr, nc, X, y, p) use norm_d, only:llf_norm implicit none integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p double precision, dimension(nc + 1), intent(out) :: r call diff_llf(llf_norm, nr, nc, X, y, p, r) end subroutine subroutine llfag(r, nr, nc, X, y, p) use norm_d, only:llfg_norm implicit none integer, intent(in) :: nr, nc double precision, dimension(nr, nc), intent(in) :: X double precision, dimension(nr), intent(in) :: y double precision, dimension(nc + 1), intent(in) :: p double precision, dimension(nc + 1), intent(out) :: r call llfg_norm(nr, nc, X, y, p, r) end subroutine
ファイル名を llf.f90 として保存します。
FortranのコードをR用のDLLにする
コマンドラインとしたら、.soか.dllができます。なお、WindowsではRtoolsが要ります。
R CMD SHLIB llf.f90 diff_llf.f90 llf_norm.f90
# RからFortranのコードを呼び出す
ベンチマーク目的
# データセットを乱数からつくる # 真のパラメーターは(1, 1, 1, -1) df01 <- with(new.env(), { n <- 500000 x <- runif(n) z <- runif(n) y <- 1 + x - z + rnorm(n) data.frame(y, x, z) }) # 推定しやすいように展開 X <- model.matrix(y ~ x + z, data = df01) nr <- nrow(X) nc <- ncol(X) p <- rep(1, nc + 1) # FortranのDLL（llf.soもしくはllf.dll）を呼び出す dll <- paste("llf", .Platform$dynlib.ext, sep="") fs <- dyn.load(dll) # 尤度関数、一階数値微分、一階解析微分 llf_f <- \(p) .Fortran("llf", double(1), nr, nc, X, df01$y, p)[[1]] llf_g <- \(p) .Fortran("llfg", double(length(p)), nr, nc, X, df01$y, p)[[1]] llf_ag <- \(p) .Fortran("llfag", double(length(p)), nr, nc, X, df01$y, p)[[1]] # 解析による微分 print(system.time({ r_ag <- optim(p, \(p) -llf_f(p), \(p) -llf_ag(p), method = "BFGS") })) # 数値微分 print(system.time({ r_ng <- optim(p, \(p) -llf_f(p), \(p) -llf_g(p), method = "BFGS") })) # グラディエントを与えない（内部的に数値微分しているはず） print(system.time({ r_wg <- optim(p, \(p) -llf_f(p), method = "BFGS") })) # DLLを開放（しないとFortranのソースコードをリコンパイルできない） dyn.unload(dll)
速度比較に使ったRのコード
lmコマンドでOLSをすれば済むので、最尤法で正規分布を推定する人はほとんどいないと思いますが、計算手順の確認および比較のために用意しました。なお、ループなし一行で書ける尤度関数はそうはないはずなので、Rのパフォーマンス評価としてはRを美化しすぎている事例になります。
# 尤度関数 R_llf <- function(p){ sum(dnorm(df01$y - X %*% p[2:(1 + nc)], sd = abs(p[1]), log = TRUE)) } # 一階微分 R_llfg <- function(p){ y <- df01$y nr <- nrow(X) nc <- ncol(X) theta <- p[1] beta <- p[-1] res <- y - X %*% beta dtheta <- -1*nr/theta + sum(res^2)/theta^3 dbeta <- numeric(nc) for(i in 1:nc){ dbeta[i] <- sum(X[, i]*2*res)/theta^2/2 } c(dtheta, dbeta) } print(system.time({ r_Ra <- optim(p, \(p) -R_llf(p), \(p) -R_llfg(p), method = "BFGS") }))
