复杂函数求导的 R 语言实践
概述
求导是求解非线性方程的重要途径之一。简单方程的求导，略懂数学知识的人手动就可以算出。然而，在实际一些算法中，求导的函数往往较为复杂。当使用如 Newton-raphson 求解时候，就需要花费大量时间来计算求导公式或导数初值。R 语言是一门统计学相关的语言，其中包含大量求解数学算法的包以及函数等，尤其方便非相关专业的编程人员的使用。R 语言开发环境下载安装简单，没有复杂的框架，运行直观容易，能够快速进行开发并运行结果。作者以个人经验为谈，旨在让初学者少走弯路，快速掌握 R 语言求导的技巧。
简单函数求导
1. Deriv 方法
deriv 是 R 中最常见的求导方法。其英文全拼“derivative”也易于使人理解。下面的例子是对公式 y=x^2 中的求导，其中须注意表达式表示方法。作为必须连接符“~”不能缺省，目标变量 y 可以省略。如不使用连接符“~”，也可采用 expression(x^2) 来直接声明表达式。deriv 方法求导后会返回 expression 类型的变量。然后通过 eval 函数赋值给 x 来计算出导数结果。举例如下：
grad <- deriv(y~x^2,"x") #deriv 方法求取一阶导数
grad #调用 grad 变量
expression({
.value <- x^2
.grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
.grad[, "x"] <- 2 * x
attr(.value, "gradient") <- .grad
.value
})
x=2 #给参数赋值
eval(grad) #给导数公式赋值
[1] 4
attr(,"gradient")
x
[1,] 4
Deriv 还可返回函数类型的结果，参数即是所求导的变量 x，当然也可以设置多个参数，这样在计算结果时可直接通过函数对参数赋值。例如：
grad <- deriv(expression(x^2),"x",func=TRUE) #func=TRUE 令返回值为 function 类型
grad(2) #grad 函数中直接参数赋值计算结果
[1] 4
attr(,"gradient")
x
[1,] 4
grad <- deriv(expression(x^2*y*z),"x",function(x,y,z){}) #用 function(){} 传递多个参数
grad(1,2,3)
[1] 6 #原函数代参结果
attr(,"gradient")
x
[1,] 12 #导函数代参结果
grad <- deriv(expression(x^2*y*z),"x",function.arg = c("x","y","z")) #用 function.arg 属性传递多个参数
grad(1,2,3)
[1] 6
attr(,"gradient")
x
[1,] 12
2. Deriv 方法
Deriv3 是求二阶导的方法，与 Deriv 使用方法相同，返回结果为二阶导数。
deriv3(y~x^4,"x")
expression({
.value <- x^4
.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"] <- 4 * x^3
.hessian[, "x", "x"] <- 4 * (3 * x^2)
attr(.value, "gradient") <- .grad
attr(.value, "hessian") <- .hessian
.value
})
3. D 方法
D 方法与 deriv 的最大不同之处在于其返回类型为“call”类型，这也意味着它更容易嵌入在其他方法中被引用。其次它能够求取二阶偏导。再者 D 方法的嵌套使用可以求取任意高阶导数。最后 D 的结果要比 deriv 的结果直观一些，deriv 常会把返回的表达式转换成复合表达式的形式，不如 D 方法显示结果直接。 D 方法的调用方式基本和 deriv 相同，在此不再赘述。二阶偏导可通过 D 方法嵌套的方式先对 x 求导再对 y 求导来获得。
#D 方法求二阶偏导
dxy = D(D(expression(x^2*y^2),"x"),"y")
dxy2 * x * (2 * y)
带参函数求导
假设带参函数如下：
其中
对上面的 y 函数进行求导，当 x 为奇数时求解第一部分，偶数时求解第二部分，所以 R 中实现方式如下。
library(Deriv)
options(digits=16)
param = c(2.0,3.0,6.0,8.0,7.0) #参数向量
dx_o = D(expression(log(x)),"x") #x 为奇数时求导公式
dx_e = D(expression(log(x^2)),"x") #x 为偶数时求导公式
dx_e
2 * x/x^2
dx_e = Simplify(dx_e) #简化公式
dx_e
2/x
grad = 0.0
for(i in 1:length(param)){ #参数循环
if(param[i]%%2){ #参数奇偶判断
grad = grad + eval(dx_e,list(x=param[i])) #参数代入相应公式进行累加求和
}else{
grad = grad + eval(dx_o,list(x=param[i]))
}
}
grad # 打印一阶导数公式计算结果
[1] 1.744047619047619
从本例中可以看出，求和在求导中最为容易，把整体函数可以拆解成单独的个体来看待，最终再合并即可。上面的 R code 引入了 eval 赋值的另外一种用法-list。当公式中包含多个参数，使用 list 方法赋值最为方便，例如 eval(expression(a+b+c)，list(a=1,b=2,c=3))。 除此之外，还使用到 Simplify 方法来简化表达式。像 dx_e 实际输出结果为 2*x/x^2，因为 D 方法或 Deriv 方法只作求导，不做简化处理，因此得到的结果可读性不强。所以 Simplify 方法对表达式进行简化处理，从而使公式更加明朗，但它并不是 R 自带方法，因此需要引入 library(Deriv) 才能使用。在引入 library 之前，需要先通过命令 install.packages(“Deriv”) 来安装该算法包。
复杂函数求导
1. 函数拼接
对上面的嵌套函数，如果全部展开来书写原函数，其一消耗时间，其二无法确保正确性，所以磨刀不误砍柴工，下面是针对公式 A 的 R code 实现。
#公式 A 嵌套拼接
epsilo <- function(p){ #拼接 epsilo
expr = "log(y)-Alpha0"
for(i in 1:p){
expr = paste(expr,"-Alpha",i,"*X",i,sep="")
}
return(expr)
}
A <- function(epsilo){ #拼接公式 A
expr = paste("exp(",epsilo,")/delta",sep="")
return(expr)
} )
p=2
D(parse(text=A(epsilo(p))),"Alpha0") #将公式 A 的字符串形式转换为表达式
-(exp(log(y) - Alpha0 - Alpha1 * X1 - Alpha2 * X2)/delta)
从上面可以看出，面对较长或嵌套的公式，我们可以简化成字符串拼接的方式，即 paste 方法来堆砌表达式，然后使用 parse(text=字符串) 将字符串转换成表达式类型，此方法在运算中非常实用。这样再复杂的函数都可以先拼接成字符串再转换成要求导的表达式，后面再做求导操作就非常容易了。
2. 函数转换
所以，公式 Y 求导可以根据上面的等式进行转换。公式 Y 转换如下：
然后，转换后的式子明显变繁为简，仅需重点求取 A关于x的偏导即可。
#函数转换求导数公式
A <- function(p){ #拼接 A 式子
alpha = c(p)
expr = "log(y) - Alpha0"
for(j in 1:p){
expr = paste(expr," - X",j,"* Alpha",j,sep = "")
}
return(expr)
}
p=2 #指定 parameter 个数
dA = deparse(D(parse(text=A(p)),"Alpha1")) #A 式对变量 Alpha1 求导，再将表达式转换成字符串
dY = NULL
dY = paste(dA,"* log(1-(",A(p),"))-(",A(p),")*1/(1-(" ,A(p),"))*", dA) #拼接 Y 的导数式子
Simplify(parse(text=dY)) #简化导数公式
expression(X1 * ((log(y) - (Alpha0 + Alpha1 * X1 + Alpha2 * X2))/(1 +
Alpha0 + Alpha1 * X1 + Alpha2 * X2 - log(y)) - log(1 + Alpha0 +
Alpha1 * X1 + Alpha2 * X2 - log(y))))
在如上代码中，deparse 用于将表达式转换成字符串再行进行拼接以完成 Y 的求导公式。所以，parse 方法与 deparse 是相对的，前者将字符串转换成表达式，后者将表达式转换成字符串。如此，复杂的表达式便可以使用子表达式转换成字符串做任意拼接。
3. 特殊函数
R 语言可以对很多特殊函数直接求导，从而使得求导更为方便。下面以 Gamma 函数与 Bessel 函数为示例说明。 (1) Gamma 函数 R 语言不仅可以直接调用 gamma 函数，而且同样可以对 gamma 函数进行求导。
#gamma 函数求导数公式
D(expression((gamma(x))^2),"x")
2 * (gamma(x) * digamma(x) * (gamma(x)))
(2) BesselK 函数 R 语言不能对 BesselK 函数直接进行求导，但是能够对 cosh 函数求导。BesselK 实际需要调用 cosh 子函数，因此，将 BesselK 转换成 cosh 相关表达式后就能使用 R 求导了。例如：
#BesselK 函数部分求导数公式
Func = expression(exp(-z*cosh(t))*cosh(v*t))
D(Func,"z")
exp(-z * cosh(t)) * (sinh(v * t) * t)
因此，BesselK(v,z) 求导后的公式为：
非线性方程求解
求导在求解非线性方程中有很普遍的应用，这里以 R 的 newtonraphson 方法中自带的例子举例。
# newtonraphson 方法求解非线性方程
ftn4 <- function(x) {
# returns function value and its derivative at x
fx <- log(x) - exp(-x)
dfx <- 1/x + exp(-x) #手动输入导数公式
return(c(fx, dfx))
}
library(spuRs)
newtonraphson(ftn4, 2, 1e-6)
At iteration 1 value of x is: 1.12202
At iteration 2 value of x is: 1.294997
At iteration 3 value of x is: 1.309709
At iteration 4 value of x is: 1.3098
Algorithm converged
[1] 1.3098
在 newtonraphson 方法中，输入包含原函数、导函数与初值。其中导函数需要从原函数求导而来。所以如采用本文中的求导方法，以上代码可更改如下：
# D 方法求导作为 newtonraphson 输入
ftn4 <- function(x) {
# returns function value and its derivative at x
fx <- expression(log(x) - exp(-x))
dfx <- D(fx,"x") #D 方法求导获取导数公式
return(c(eval(fx), eval(dfx)))
}
newtonraphson(ftn4, 2, 1e-6)
At iteration 1 value of x is: 1.12202
At iteration 2 value of x is: 1.294997
At iteration 3 value of x is: 1.309709
At iteration 4 value of x is: 1.3098
Algorithm converged
[1] 1.3098
本例中，使用 D 方法求导来代替手动输入导数公式。这样，仅保证原函数输入正确即可，导函数使用 D 方法或 Deriv 方法求出，从而省去计算导函数并验证导函数是否正确的步骤。在实际应用过程中，原函数已经相当复杂，如若计算导函数并键入公式，耗时费力且结果未必正确，因此使用 R 的求导方法可以大大简化该繁琐步骤。
结束语
求导在数学算法应用中相当普遍，例如线性回归模型，Cox 回归模型，参数化回归模型中等等。本文旨在利用 R 中的求导方式使得该流程简单化，省去费力耗时的工作，从而达到事半功倍的效果。本文是在作者实际使用摸索的基础上摘取一些可以帮助初学者的例子，希望可以使需要的人受益。
欢迎大家关注微信公众号：数萃大数据
深度学习培训班【上海站】
时间：2017年12月23-24日
地点：上海创梦云实训创新中心
更多详情，请扫描下面二维码
相关内容
微信扫码咨询专知VIP会员