R语言编程与仿真
作为此书的笔记, 全部代码在此
挑选比较值得记录的地方作为一个笔记吧
最优化
牛顿法
这是最简单的, 最为牛顿法寻根的一个扩展, 只需要把寻根迭代公式运用到导数求根的过程中就可以了, 相比最速下降法的优势是有二阶导数信息,因此收敛速度会更快一些, 但也有不收敛的可能$$x(n+1)=x(n)-\frac{f’(x(n))}{f’’(x(n))}$$
newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
# Newton's method for optimisation, starting at x0
# f3 is a function that given x returns the vector
# (f(x), f'(x), f''(x)), for some f
x <- x0
f3.x <- f3(x)
n <- 0
while ((abs(f3.x[2]) > tol) & (n < n.max)) {
x <- x - f3.x[2]/f3.x[3]
f3.x <- f3(x)
n <- n + 1
}
if (n == n.max) {
cat('newton failed to converge\n')
} else {
return(x)
}
}
gamma.2.3 <- function(x) {
# gamma(2,3) density
if (x < 0) return(c(0, 0, 0))
if (x == 0) return(c(0, 0, NaN))
y <- exp(-2*x)
return(c(4*x^2*y, 8*x*(1-x)*y, 8*(1-2*x^2)*y))
}
黄金分割法
此方法不需要一阶导数, 类似于夹逼法, 需要三个起始点进行迭代.
gsection <- function(ftn, x.l, x.r, x.m, tol = 1e-9) {
# applies the golden-section algorithm to maximise ftn
# we assume that ftn is a function of a single variable
# and that x.l < x.m < x.r and ftn(x.l), ftn(x.r) <= ftn(x.m)
#
# the algorithm iteratively refines x.l, x.r, and x.m and terminates
# when x.r - x.l <= tol, then returns x.m
# golden ratio plus one
gr1 <- 1 + (1 + sqrt(5))/2
# successively refine x.l, x.r, and x.m
f.l <- ftn(x.l)
f.r <- ftn(x.r)
f.m <- ftn(x.m)
while ((x.r - x.l) > tol) {
if ((x.r - x.m) > (x.m - x.l)) {
y <- x.m + (x.r - x.m)/gr1
f.y <- ftn(y)
if (f.y >= f.m) {
x.l <- x.m
f.l <- f.m
x.m <- y
f.m <- f.y
} else {
x.r <- y
f.r <- f.y
}
} else {
y <- x.m - (x.m - x.l)/gr1
f.y <- ftn(y)
if (f.y >= f.m) {
x.r <- x.m
f.r <- f.m
x.m <- y
f.m <- f.y
} else {
x.l <- y
f.l <- f.y
}
}
}
return(x.m)
}
多元优化问题
考虑多元优化问题, 需要引进梯度和黑塞矩阵两个因素
局部最小极值的必要条件是$${\nabla}f=0$$
充分条件是$${\nabla}f(x)=0,且H(x)是正定的$$所正定, f(x)沿任何方向的曲率都是正的(对于局部最大值有相反的结论, 即负定.
最速上升法
注: 这本书讨论的都是求极大值, 不过我感觉最小值习惯一点, 所以把他讨论的内容改成了极小值的条件.
迭代公式是$$x(n+1)=x(n)-{\alpha}{\nabla}f(x(n))$$, 其中$\alpha$是正数, 最为一个参数, 这个问题可以化为一个单参数优化问题, 找到最优的$\alpha$使得下一次迭代的函数值最小.即$$g(\alpha)=f(x(n)-{\alpha}{\nabla}f(x(n)))$$
source("../scripts/linesearch.r")
ascent <- function(f, grad.f, x0, tol = 1e-9, n.max = 100) {
# steepest ascent algorithm
# find a local max of f starting at x0
# function grad.f is the gradient of f
x.old <- x0
x <- line.search(f, x0, grad.f(x0))
n <- 1
while ((f(x) - f(x.old) > tol) & (n < n.max)) {
x.old <- x
x <- line.search(f, x, grad.f(x))
n <- n + 1
}
return(x)
}
线性搜索
注意这里的linesearch函数, 其用处就是找出最优的$\alpha$, 在这里使用黄金分割进行寻找
line.search <- function(f, x, y, tol = 1e-9, a.max = 2^5) {
# f is a real function that takes a vector of length d
# x and y are vectors of length d
# line.search uses gsection to find a >= 0 such that
# g(a) = f(x + a*y) has a local maximum at a,
# within a tolerance of tol
# if no local max is found then we use 0 or a.max for a
# the value returned is x + a*y
if (sum(abs(y)) == 0) return(x) # g(a) constant
g <- function(a) return(f(x + a*y))
# find a triple a.l < a.m < a.r such that
# g(a.l) <= g(a.m) and g(a.m) >= g(a.r)
# a.l
a.l <- 0
g.l <- g(a.l)
# a.m
a.m <- 1
g.m <- g(a.m)
while ((g.m < g.l) & (a.m > tol)) {
a.m <- a.m/2
g.m <- g(a.m)
}
# if a suitable a.m was not found then use 0 for a
if ((a.m <= tol) & (g.m < g.l)) return(x)
# a.r
a.r <- 2*a.m
g.r <- g(a.r)
while ((g.m < g.r) & (a.r < a.max)) {
a.m <- a.r
g.m <- g.r
a.r <- 2*a.m
g.r <- g(a.r)
}
# if a suitable a.r was not found then use a.max for a
if ((a.r >= a.max) & (g.m < g.r)) return(x + a.max*y)
# apply golden-section algorithm to g to find a
a <- gsection(g, a.l, a.r, a.m)
return(x + a*y)
}
由于黄金分割需要三个起始点, 所以这里首先进行搜寻起始点, 当搜索到后, 调用之前的黄金分割找出$\alpha$, 计算出下一步的$x(n+1)$, 返回给最速下降函数进行迭代.
多维下的牛顿法
最速上升法使用了梯度信息, 但是如果使用二阶信息黑丝矩阵, 可以构造出更快的方法, 这里就引出了牛顿法, 对于比较接近的两个点, $y$与$x$我们有$$f(y)=f(x)+(y-x)^{T}{\nabla}f(x)+\frac{1}{2}(y-x)^{T}H(x)(y-x)$$, 即泰勒展开的二阶近似.对两边求导$${\nabla}f(y)={\nabla}f(x)+H(x)(y-x)$$, 令梯度为0, 我们就有了迭代公式$$x(n+1)=x(n)-H(x(n))^{-1}{\nabla}f(x(n))$$
这里如果H(x)是奇异的(没有逆)那么就无法继续迭代, 并且就算不是奇异的也不乏保证收敛性, 这里需要满足的条件是$H(x)$是正定矩阵, 且满足$h(x)$利普希茨连续.
有如下代码
这里输入就是一个可以计算x处函数值, 梯度, 黑塞矩阵的函数, 初始值, (误差最小, 迭代限制)
newton <- function(f3, x0, tol = 1e-9, n.max = 100) {
# Newton's method for optimisation, starting at x0
# f3 is a function that given x returns the list
# {f(x), grad f(x), Hessian f(x)}, for some f
x <- x0
f3.x <- f3(x)
n <- 0
while ((max(abs(f3.x[[2]])) > tol) & (n < n.max)) {
x <- x - solve(f3.x[[3]], f3.x[[2]])
f3.x <- f3(x)
n <- n + 1
}
if (n == n.max) {
cat('newton failed to converge\n')
} else {
return(x)
}
}
R中有可以用来求梯度和黑塞矩阵的函数deriv
如
df <- deriv(z ~ sin(x^2/2-y^2/4)*cos(2*x-exp(y)),c('x','y'),func=TRUE,hessian=TRUE)
f3 <- function(x){
dfx <- df(x[1],x[2])
f <- dfx[1]
gradf <- attr(dfx,'gradient')[1,]
hessianf <- attr(dfx,'hessian')[1,,]
return(list(f,gradf,hessianf))
}
R内部的函数
一般的有optimize和optim,其中optimize(黄金分割法和抛物线插值法的结合)对应的一维度情况, optim对应的多维情况(使用的是NM, 伪牛顿法, 共轭梯度法), 一般不适用牛顿法和最速下降法(主要是计算速度太慢了)
对于求根
方法有
- 牛顿法
- 割线法(牛顿法的近似, 不用求导数)
- 二分法(不用求导数)