第 24 章 数值优化

数值微分、,带约束与无约束,多元函数优化,带约束和无约束,整数优化 凸优化 Anqi Fu 开发的 CVXR Jelmer Ypma 开发的 nloptr 和 Berwin A. Turlach 开发的 quadprog

24.1 一元非线性优化

复合函数求极值

\[ g(x) = \int_{0}^{x} -\sqrt{t}\exp(-t^2) dt, \\ f(y) = \int_{0}^{y} g(s) \exp(-s) ds \]

g <- function(x) {
  integrate(function(t) {
    -sqrt(t) * exp(-t^2)
  }, lower = 0, upper = x)$value
}

f <- function(y) {
  integrate(function(s) {
    Vectorize(g, "x")(s) * exp(-s)
  }, lower = 0, upper = y)$value
}

optimize(f, interval = c(10, 100), maximum = FALSE)
## $minimum
## [1] 66.84459
## 
## $objective
## [1] -0.3201572

计算积分的时候,输入了一系列 s 值,参数是向量,而函数 g 只支持输入参数是单个值,g(c(1,2)) 则会报错

g(1)
## [1] -0.453392

24.2 多元无约束非线性优化

GA 实现遗传算法

24.3 凸二次规划

在 R 中使用 quadprog 包求解二次规划24,而 ipoptr 包可用来求解一般的非线性约束的非线性规划25,quadprogXT 包用来求解带绝对值约束的二次规划,pracma 包提供 quadprog 函数就是对 quadprog 包的 solve.QP 进行封装,使得调用风格更像 Matlab 而已。quadprog 包实现了 Goldfarb and Idnani (1982, 1983) 提出的对偶方法,主要用来求解带线性约束的严格凸二次规划问题。

\[\min(-d^{T}b+1/2b^{T}Db), A^{T}b \geq b_{0}\]

solve.QP(Dmat, dvec, Amat, bvec, meq = 0, factorized = FALSE)

参数 DmatdvecAmatbvec 分别对应二次规划问题中的 \(D,d,A,b_{0}\)

现在有如下二次规划问题

\[ D = 2* \begin{bmatrix}1 & -1/2\\ -1/2 & 1 \end{bmatrix}, d = (-3,2), A = \begin{bmatrix}1 & 1\\ -1 & 1 \\ 0 & -1 \end{bmatrix}, b_{0} = (2,-2,-3) \]

上述二次规划问题的可行域如图所示

plot(0, 0,
  xlim = c(-2, 5.5), ylim = c(-1, 3.5), type = "n",
  xlab = "x", ylab = "y", main = "Feasible Region"
)
polygon(c(2, 5, -1), c(0, 3, 3), border = TRUE, lwd = 2, col = "gray")
可行域

图 24.1: 可行域

quadprog 包的 solve.QP 函数求解二次规划

library(quadprog)
Dmat <- 2 * matrix(c(1, -1 / 2, -1 / 2, 1), nrow = 2, byrow = TRUE)
dvec <- c(-3, 2)
A <- matrix(c(1, 1, -1, 1, 0, -1), ncol = 2, byrow = TRUE)
bvec <- c(2, -2, -3)
Amat <- t(A)
sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = 0)
sol
## $solution
## [1] 0.1666667 1.8333333
## 
## $value
## [1] -0.08333333
## 
## $unconstrained.solution
## [1] -1.3333333  0.3333333
## 
## $iterations
## [1] 2 0
## 
## $Lagrangian
## [1] 1.5 0.0 0.0
## 
## $iact
## [1] 1

在可行域上画出等高线,表示目标解的位置,图中红点表示无约束下的解,黄点表示线性约束下的解

qp_sol <- sol$solution # 二次规划的解
uc_sol <- sol$unconstrained.solution # 无约束情况下的解
# 画图
library(lattice)
x <- seq(-2, 5.5, length.out = 500)
y <- seq(-1, 3.5, length.out = 500)
grid <- expand.grid(x = x, y = y)
grid$z <- with(grid, x^2 + y^2 - x * y + 3 * x - 2 * y + 4)
levelplot(z ~ x * y, grid,
  cuts = 40,
  panel = function(...) {
    panel.levelplot(...)
    panel.polygon(c(2, 5, -1), c(0, 3, 3),
      border = TRUE,
      lwd = 2, col = "transparent"
    )
    panel.points(c(uc_sol[1], qp_sol[1]),
      c(uc_sol[2], qp_sol[2]),
      lwd = 5, col = c("red", "yellow"), pch = 19
    )
  },
  colorkey = TRUE,
  col.regions = terrain.colors(40)
)
无约束和有约束条件下的解

图 24.2: 无约束和有约束条件下的解

可行域 polypath 线性约束 非线性约束如何

24.4 对数似然

set.seed(1234)
n <- 20 # 随机数的个数
x <- rexp(n, rate = 5) # 服从指数分布的随机数
m <- 40 # 网格数
mu <- seq(mean(x) - 1.5 * sd(x) / sqrt(n),
          mean(x) + 1.5 * sd(x) / sqrt(n),
          length.out = m
)
sigma <- seq(0.8 * sd(x), 1.5 * sd(x), length.out = m)
tmp <- expand.grid(x = mu, y = sigma)
loglikelihood <- function(b) -sum(dnorm(x, b[1], b[2], log = TRUE))
pp <- apply(tmp, 1, loglikelihood)
z <- matrix(pp, m, m)
nbcol <- 100
color <- hcl.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -m] + z[-m, -1] + z[-m, -m]
facetcol <- cut(zfacet, nbcol)

# cairo_pdf(file = "log-lik.pdf",width = 6,height = 5)
par(mar = c(0.1, 2, 0.1, 0.1))
persp(mu, sigma, z,
      xlab = "\n \u03bc", ylab = "\n \u03c3",
      zlab = "\n log-likelihood",
      border = NA,
      ticktype = "simple",
      col = color[facetcol],
      theta = 50, phi = 25,
      r = 60, d = 0.1, expand = .6,
      ltheta = 90, lphi = 180,
      shade = 0.1, nticks = 5 # , box = TRUE,axes = TRUE
)

# 添加极大值点
# 除指数分布外,还有正态、二项、泊松分布观察其似然曲面的特点,都是单峰,有唯一极值点
# 再考虑正态混合模型的似然曲面
# dev.off()
# ```{r log-lik,fig.cap="对数似然函数曲面",fig.asp=1,out.width="70%",dev="png",dev.args=list(type="cairo", bg = "transparent"), dpi = 300}
# 
# ```
# Sys.setenv(R_GSCMD = "C:/Program Files/gs/gs9.26/bin/gswin64c.exe")
# embedFonts(file = 'log-lik.pdf',outfile = 'log-lik-emfont.pdf',
#            fontpaths = system.file("fonts", package = "fontcm"))

24.5 一元非线性方程

牛顿-拉弗森方法

24.6 微分方程

ode45 求解偏微分方程

library(nonlinearTseries)
library(plot3D)
lor <- lorenz(do.plot = F)

scatter3D(lor$x, lor$y, lor$z,
  ann = FALSE, col = terrain.colors(25),
  type = "o", cex = 0.3,
  colkey = FALSE, box = FALSE
)