R optim 通用优化

说明

基于Nelder-Mead、quasi-Newton和共轭梯度算法的通用优化。它包括用于盒约束优化和模拟退火的选项。

用法

optim(par, fn, gr = NULL, ...,
      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",
                 "Brent"),
      lower = -Inf, upper = Inf,
      control = list(), hessian = FALSE)

optimHess(par, fn, gr = NULL, ..., control = list())

参数

par	要优化的参数的初始值。
fn	要最小化(或最大化)的函数，第一个参数是要进行最小化的参数向量。它应该返回一个标量结果。
gr	返回 "BFGS" 、 "CG" 和 "L-BFGS-B" 方法的梯度的函数。如果是NULL，将使用有限差分近似。 对于"SANN" 方法，它指定一个函数来生成新的候选点。如果是NULL，则使用默认的高斯马尔可夫内核。
...	要传递给 fn 和 gr 的更多参数。
method	要使用的方法。查看具体信息'。可以缩写。
lower, upper	"L-BFGS-B" 方法的变量界限，或搜索方法 "Brent" 的界限。
control	控制参数list。查看具体信息'。
hessian	逻辑性强。是否应该返回数值微分的 Hessian 矩阵？

细节

请注意，... 之后的参数必须完全匹配。

默认情况下，optim 执行最小化，但如果 control$fnscale 为负数，则会最大化。 optimHess 是一个辅助函数，用于在忘记hessian = TRUE 的情况下在后期计算 Hessian 矩阵。

默认方法是 Nelder 和 Mead (1965) 方法的实现，它仅使用函数值，并且很健壮，但相对较慢。它对于不可微函数来说效果相当好。

方法 "BFGS" 是 quasi-Newton 方法(也称为可变度量算法)，具体由 Broyden、Fletcher、Goldfarb 和 Shanno 于 1970 年同时发布。它使用函数值和梯度来构建要优化的表面的图片。

方法 "CG" 是基于 Fletcher 和 Reeves (1964) 的共轭梯度方法(但可以选择 Polak-Ribiere 或 Beale-Sorenson 更新)。共轭梯度方法通常比 BFGS 方法更脆弱，但由于它们不存储矩阵，因此它们可能在更大的优化问题中取得成功。

方法"L-BFGS-B" 是 Byrd 等人的方法。等人。 (1995) 允许框约束，即每个变量都可以给出一个下限和/或上限。初始值必须满足约束条件。这使用 BFGS quasi-Newton 方法的有限内存修改。如果提供了非平凡边界，则将选择此方法，并带有警告。

Nocedal and Wright (1999)是对前三种方法的综合参考。

默认情况下，方法 "SANN" 是 Belisle (1992) 中给出的模拟退火的变体。模拟退火属于随机全局优化方法类。它仅使用函数值，但速度相对较慢。它也适用于不可微函数。此实现使用 Metropolis 函数来计算接受概率。默认情况下，下一个候选点是从高斯马尔可夫核生成的，其比例与实际温度成比例。如果给出生成新候选点的函数，方法"SANN"也可以用于解决组合优化问题。温度根据 Belisle (1992, p. 890) 中给出的对数冷却时间表降低；具体来说，温度设置为 temp / log(((t-1) %/% tmax)*tmax + exp(1)) ，其中 t 是当前迭代步骤，temp 和 tmax 可通过 control 指定，请参见下文。请注意，"SANN" 方法很大程度上取决于控制参数的设置。它不是通用方法，但对于在非常粗糙的表面上获得良好的值非常有用。

方法 "Brent" 仅适用于一维问题，使用 optimize() 。当 optim() 在只能指定 method 的其他函数中使用时，它会很有用，例如在 stats4 包中的 mle 中。

如果函数 fn 无法按提供的值进行计算，则该函数可以返回 NA 或 Inf ，但初始值必须具有可计算的有限值 fn 。 (除了方法 "L-BFGS-B" 之外，该方法的值应始终是有限的。)

optim 可以递归使用，并且可以用于单个参数，也可以用于多个参数。它还接受零长度 par ，并仅使用该参数计算函数。

control 参数是一个可以提供以下任何组件的列表：

trace

非负整数。如果是肯定的，则会生成有关优化进度的跟踪信息。较高的值可能会产生更多的跟踪信息：对于方法"L-BFGS-B"，有六个跟踪级别。 (要准确理解这些内容，请参阅源代码：更高级别提供更多详细信息。)

fnscale

优化期间应用于 fn 和 gr 值的整体缩放。如果为负，则将问题转化为最大化问题。对 fn(par)/fnscale 进行优化。

parscale

参数缩放值的向量。优化是在 par/parscale 上执行的，从任何元素中的单位变化都会产生缩放值的单位变化的意义上来说，这些应该具有可比性。 method = "Brent" 不使用(也不需要)。

ndeps

梯度有限差分近似的步长向量，采用 par/parscale 尺度。默认为 1e-3 。

maxit

最大迭代次数。对于基于导数的方法，默认为 100；对于 "Nelder-Mead"，默认为 500。

对于 "SANN" maxit 给出函数评估的总数：没有其他停止标准。默认为 10000 。

abstol

绝对收敛容差。仅对非负函数有用，作为达到零的容差。

reltol

相对收敛容差。如果无法将值一步减少 reltol * (abs(val) + reltol) 倍，算法就会停止。默认为 sqrt(.Machine$double.eps) ，通常约为 1e-8 。

alpha , beta , gamma

"Nelder-Mead" 方法的缩放参数。 alpha 是反射系数(默认 1.0)，beta 是收缩系数(0.5)，gamma 是膨胀系数(2.0)。

REPORT

如果 control$trace 为正，则 "BFGS" 、 "L-BFGS-B" 和 "SANN" 方法的报告频率。对于 "BFGS" 和 "L-BFGS-B" 默认为每 10 次迭代，对于 "SANN" 默认为每 100 个温度。

warn.1d.NelderMead

logical 指示(默认)"Nelder-Mead" 方法在用于一维最小化时是否应发出警告信号。由于警告有时是不适当的，您可以通过将此选项设置为 false 来抑制它。

type

对于共轭梯度法。对于 Fletcher-Reeves 更新，采用值 1；对于 Polak-Ribiere，采用值 2；对于 Beale-Sorenson，采用值 3。

lmm

是一个整数，给出 "L-BFGS-B" 方法中保留的 BFGS 更新数量，默认为 5 。

factr

控制"L-BFGS-B"方法的收敛。当目标的减少量在机器公差的这个因子内时，就会发生收敛。默认值为 1e7 ，即大约 1e-8 的容差。

pgtol

有助于控制 "L-BFGS-B" 方法的收敛。它是当前搜索方向上投影梯度的容差。当检查被抑制时，该值默认为零。

temp

控制"SANN"方法。它是冷却计划的起始温度。默认为 10 。

tmax

是 "SANN" 方法在每个温度下的函数评估次数。默认为 10 。

赋予 par 的任何名称都将复制到传递给 fn 和 gr 的向量。请注意，par 的其他属性不会被复制。

传递给 fn 的参数向量具有特殊语义，并且可以在调用之间共享：函数不应更改或复制它。

值

对于 optim ，包含组件的列表：

par	找到的最佳参数集。
value	fn 对应于 par 的值。
counts	一个二元素整数向量，分别给出对 fn 和 gr 的调用次数。这不包括计算 Hessian 矩阵(如果需要)所需的调用，以及对 fn 计算梯度有限差分近似值的任何调用。
convergence	整数代码。 0 表示成功完成("SANN" 和 "Brent" 始终是这种情况)。可能的错误代码有 1 表示已达到迭代限制maxit。 10 表示 Nelder-Mead 单纯形的简并性。 51 表示来自 "L-BFGS-B" 方法的警告；有关更多详细信息，请参阅组件message。 52 指示 "L-BFGS-B" 方法中的错误；有关更多详细信息，请参阅组件message。
message	给出优化器返回的任何附加信息的字符串，或 NULL 。
hessian	仅当参数 hessian 为 true 时。对称矩阵给出所找到解的 Hessian 矩阵估计。请注意，即使框约束处于活动状态，这也是无约束问题的 Hessian 矩阵。

对于 optimHess ，适用 hessian 组件的说明。

注意

optim 将与一维 par 一起使用，但默认方法不能很好地工作(并且会发出警告)。方法 "Brent" 使用 optimize 并且需要边界可用；如果不是的话，"BFGS" 通常也能正常工作。

例子

require(graphics)

fr <- function(x) {   ## Rosenbrock Banana function
    x1 <- x[1]
    x2 <- x[2]
    100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}
grr <- function(x) { ## Gradient of 'fr'
    x1 <- x[1]
    x2 <- x[2]
    c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
       200 *      (x2 - x1 * x1))
}
optim(c(-1.2,1), fr)
(res <- optim(c(-1.2,1), fr, grr, method = "BFGS"))
optimHess(res$par, fr, grr)
optim(c(-1.2,1), fr, NULL, method = "BFGS", hessian = TRUE)
## These do not converge in the default number of steps
optim(c(-1.2,1), fr, grr, method = "CG")
optim(c(-1.2,1), fr, grr, method = "CG", control = list(type = 2))
optim(c(-1.2,1), fr, grr, method = "L-BFGS-B")

flb <- function(x)
    { p <- length(x); sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) }
## 25-dimensional box constrained
optim(rep(3, 25), flb, NULL, method = "L-BFGS-B",
      lower = rep(2, 25), upper = rep(4, 25)) # par[24] is *not* at boundary


## "wild" function , global minimum at about -15.81515
fw <- function (x)
    10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80
plot(fw, -50, 50, n = 1000, main = "optim() minimising 'wild function'")

res <- optim(50, fw, method = "SANN",
             control = list(maxit = 20000, temp = 20, parscale = 20))
res
## Now improve locally {typically only by a small bit}:
(r2 <- optim(res$par, fw, method = "BFGS"))
points(r2$par,  r2$value,  pch = 8, col = "red", cex = 2)

## Combinatorial optimization: Traveling salesman problem
library(stats) # normally loaded

eurodistmat <- as.matrix(eurodist)

distance <- function(sq) {  # Target function
    sq2 <- embed(sq, 2)
    sum(eurodistmat[cbind(sq2[,2], sq2[,1])])
}

genseq <- function(sq) {  # Generate new candidate sequence
    idx <- seq(2, NROW(eurodistmat)-1)
    changepoints <- sample(idx, size = 2, replace = FALSE)
    tmp <- sq[changepoints[1]]
    sq[changepoints[1]] <- sq[changepoints[2]]
    sq[changepoints[2]] <- tmp
    sq
}

sq <- c(1:nrow(eurodistmat), 1)  # Initial sequence: alphabetic
distance(sq)
# rotate for conventional orientation
loc <- -cmdscale(eurodist, add = TRUE)$points
x <- loc[,1]; y <- loc[,2]
s <- seq_len(nrow(eurodistmat))
tspinit <- loc[sq,]

plot(x, y, type = "n", asp = 1, xlab = "", ylab = "",
     main = "initial solution of traveling salesman problem", axes = FALSE)
arrows(tspinit[s,1], tspinit[s,2], tspinit[s+1,1], tspinit[s+1,2],
       angle = 10, col = "green")
text(x, y, labels(eurodist), cex = 0.8)

set.seed(123) # chosen to get a good soln relatively quickly
res <- optim(sq, distance, genseq, method = "SANN",
             control = list(maxit = 30000, temp = 2000, trace = TRUE,
                            REPORT = 500))
res  # Near optimum distance around 12842

tspres <- loc[res$par,]
plot(x, y, type = "n", asp = 1, xlab = "", ylab = "",
     main = "optim() 'solving' traveling salesman problem", axes = FALSE)
arrows(tspres[s,1], tspres[s,2], tspres[s+1,1], tspres[s+1,2],
       angle = 10, col = "red")
text(x, y, labels(eurodist), cex = 0.8)

## 1-D minimization: "Brent" or optimize() being preferred.. but NM may be ok and "unavoidable",
## ----------------   so we can suppress the check+warning :
system.time(rO <- optimize(function(x) (x-pi)^2, c(0, 10)))
system.time(ro <- optim(1, function(x) (x-pi)^2, control=list(warn.1d.NelderMead = FALSE)))
rO$minimum - pi # 0 (perfect), on one platform
ro$par - pi     # ~= 1.9e-4    on one platform
utils::str(ro)

来源

方法 "Nelder-Mead" 、 "BFGS" 和 "CG" 的代码最初基于 Nash (1990) 中的 Pascal 代码，由 p2c 翻译，然后进行手工优化。纳什博士同意可以免费提供该代码。

方法的代码"L-BFGS-B"基于 Zhu、Byrd、Lu-Chen 和从 Netlib 获得的 Nocedal 的 Fortran 代码(文件‘选择/lbfgs_bcm.shar’：另一个版本在‘汤姆斯/778’)。

方法 "SANN" 的代码由 A. Trapletti 贡献。

参考

Belisle, C. J. P. (1992). Convergence theorems for a class of simulated annealing algorithms on R^d. Journal of Applied Probability, 29, 885-895. doi:10.2307/3214721.

Byrd, R. H., Lu, P., Nocedal, J. and Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16, 1190-1208. doi:10.1137/0916069.

Fletcher, R. and Reeves, C. M. (1964). Function minimization by conjugate gradients. Computer Journal 7, 148-154. doi:10.1093/comjnl/7.2.149.

Nash, J. C. (1990). Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation. Adam Hilger.

Nelder, J. A. and Mead, R. (1965). A simplex algorithm for function minimization. Computer Journal, 7, 308-313. doi:10.1093/comjnl/7.4.308.

Nocedal, J. and Wright, S. J. (1999). Numerical Optimization. Springer.

也可以看看

nlm、nlminb。

optimize 用于一维最小化，constrOptim 用于约束优化。