最近筆者有計算一個函數的Gradient和Hessian的需求，所以研究了一下R 中進行這類計算的功能和套件。
以下我會介紹R 中進行數值微分和代數微分的功能。
我第一個介紹的是數值微分套件：numDeriv。他的用法很簡單：
gradgrad函數可以拿來計算函數在指定點的gradient：
library(numDeriv)
f <- function(x) x^2
grad(f, 2)
## [1] 4
grad也可以計算向量函數的gradient：
f <- function(x) x[1] * x[2]
grad(f, c(1, 2))
## [1] 2 1
hessian可以計算函數在指定點的hessian
hessian(f, 1:2)
## [,1] [,2]
## [1,] 1.684e-13 1.00e+00
## [2,] 1.000e+00 4.21e-14
由於是數值解的關係，所以可以看到二階微分會算出很小的值。
deriv接下來我要介紹一個做代數微分的功能。
首先我們把剛剛的f <- function(x) x^2寫成expression：
f <- function(x) x^2
f.exp <- expression(x^2)
接著使用內建的deriv：
deriv(f.exp, namevec = "x")
## expression({
## .value <- x^2
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
## .grad[, "x"] <- 2 * x
## attr(.value, "gradient") <- .grad
## .value
## })
R 就會回傳出計算gradient的方式。將它複製貼上建立函數，並且修改內容讓它是vectorized：
g <- function(x) {
.value <- x^2
.grad <- 2 * x
.grad
}
和numDeriv的結果做個比較：
all.equal(g(2), grad(f, 2))
## [1] TRUE
deriv的另一個參數hessian設為TRUE的時候，回傳值會附加hessian的算式。
根據Bowling et al. (2009)，以下函數是standard normal cdf的logistic approximation：
\[ F(x) = \frac{1}{1+e^{-0.07056 x^3 - 1.5976 x}} \]
F.0 <- function(x) 1/(1 + exp(-0.07056 * x^3 - 1.5976 * x))
curve(pnorm, -3, 3, lty = 1)
curve(F.0, add = TRUE, col = 2, lty = 2)
肉眼根本分不出來差異！
如果我們用它來取代某些Likelihood中的standard normal cdf(當資料有censoring的時候)，在計算MLE的時候就會需要計算\( F'(x) \)和\( F''(x) \)。他們分別是：
好醜好難寫。
這時候就可以利用deriv來建立gradient and hessian。先建立gradient：
F.exp <- expression(1/(1 + exp(-0.07056 * x^3 - 1.5976 * x)))
deriv(F.exp, namevec = "x")
## expression({
## .expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
## .expr7 <- 1 + .expr6
## .value <- 1/.expr7
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
## .grad[, "x"] <- .expr6 * (0.07056 * (3 * x^2) + 1.5976)/.expr7^2
## attr(.value, "gradient") <- .grad
## .value
## })
稍作修改就可以得到：
F.1 <- function(x) {
.expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
.expr7 <- 1 + .expr6
.value <- 1/.expr7
.expr6 * (0.07056 * (3 * x^2) + 1.5976)/.expr7^2
}
再參考deriv的hessian算式：
deriv(F.exp, namevec = "x", hessian = TRUE)
## expression({
## .expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
## .expr7 <- 1 + .expr6
## .expr12 <- 0.07056 * (3 * x^2) + 1.5976
## .expr13 <- .expr6 * .expr12
## .expr14 <- .expr7^2
## .value <- 1/.expr7
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
## .hessian <- array(0, c(length(.value), 1L, 1L), list(NULL,
## c("x"), c("x")))
## .grad[, "x"] <- .expr13/.expr14
## .hessian[, "x", "x"] <- (.expr6 * (0.07056 * (3 * (2 * x))) -
## .expr13 * .expr12)/.expr14 + .expr13 * (2 * (.expr13 *
## .expr7))/.expr14^2
## attr(.value, "gradient") <- .grad
## attr(.value, "hessian") <- .hessian
## .value
## })
可以寫出：
F.2 <- function(x) {
.expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
.expr7 <- 1 + .expr6
.expr12 <- 0.07056 * (3 * x^2) + 1.5976
.expr13 <- .expr6 * .expr12
.expr14 <- .expr7^2
1/.expr7
(.expr6 * (0.07056 * (3 * (2 * x))) - .expr13 * .expr12)/.expr14 + .expr13 *
(2 * (.expr13 * .expr7))/.expr14^2
}
和numDeriv比較一下：
x <- rnorm(1)
grad(F.0, x)
## [1] 0.04958
F.1(x)
## [1] 0.04958
hessian(F.0, x)
## [,1]
## [1,] -0.1008
F.2(x)
## [1] -0.1008
完全一致！
但是用deriv建立的函數效能絕對比較好：
library(microbenchmark)
microbenchmark(grad(F.0, x), F.1(x))
## Unit: microseconds
## expr min lq median uq max neval
## grad(F.0, x) 91.989 161.198 163.853 169.821 271.03 100
## F.1(x) 4.322 7.819 8.266 9.645 38.97 100
microbenchmark(hessian(F.0, x), F.2(x))
## Unit: microseconds
## expr min lq median uq max neval
## hessian(F.0, x) 203.193 207.84 215.511 244.846 339.99 100
## F.2(x) 6.878 7.97 9.348 9.898 16.66 100
以上內容給大家參考囉！