【推荐系列】R 语言中的最优化方法

作者 | immanuel

来源 | 金融投资智库

什么是最优化方法

最优化方法与运筹学

最优化方法是在所有的可行方案中找出最优方案的方法；也称为运筹学方法。

最优化方法与数据分析

最优化方法与统计学、机器学习等数据分析方法不同，并不是从历史数据中发现规律和建立模型；最优化方法与数据分析方法的相似之处在于流程，都需要从实际问题中抽象出数学模型，然后根据数据和问题来建模，并使用计算机求解；最优化方法是数据科学中不可或缺的关键一环，很多统计模型和机器学习算法的理解与实现都需要仰仗最优化方法。

线性规划

lpSolve 是R 中最常用的线性规划工具，基于C 语言；

Rglpk、Rmosek、Rcplex 等是其他工具的R 语言接口。

非线性规划

无约束的非线性规划问题比较简单，R 基础包中包含一些很

方便的函数；

Rsolnp 和alabama 用来求解约束非线性规划的问题；

对于二次规划问题，有quadprog 等专门的R 包可以解决。

整数规划

简单的整数规划问题可以使用线性规划工具解决；

混合线性整数规划可以使用Rsymphony 包来解决，该包基于

COIN-OR 上的SYMPHONY 项目；

非线性混合整数规划问题暂无方便的R 包可以直接使用。

图优化

igraph 是sna 常用的图与网络优化的工具包；

这也是进来流行的社交网络分析的基础方法。

一个标准的最优化模型

两变量的问题可以直观地在直角坐标性中理解:

lpSolve 包的lp 函数

direction：优化的方向，默认是“min”，也可以是“max”

objective.in：目标函数的系数，向量的长度必须等于变量的

个数

const.mat：约束条件的系数矩阵

const.dir：约束的方向

const.rhs：不等式右边的数值向量

Rglpk 包

AMPL 的子集

Rglpk_read_file 可以读入模型

Rglpk_solve_LP 可以求解模型

代码示例：

library(lpSolve)

f.obj <- c(-2, -1)

f.con <- matrix (c(-3, -4, -1, 2, 1, 0, 0, 1), nrow = 4,

byrow = TRUE)

f.dir <- c(">=", ">=", ">=", ">=")

f.rhs <- c(-12, -2, 0, 0)

res <- lp("min", f.obj, f.con, f.dir, f.rhs)

res

## Success: the objective function is -7

class(res)

## [1] "lp"

res$solution

## [1] 3.2 0.6

IP

纯粹的整数规划

所有变量都必须取值为整数

所有变量都为0-1 时称为0-1 规划

MILP

混合整数线性规划

部分变量为整数

MINLP

混合整数规划

目标函数和约束条件可能同时包含整数项和非线性项

Rglpk 包示例

library(Rglpk)

mod.file <- system.file("examples", "optimization",

"sudoku.mod", package = "rinds")

mod <- Rglpk_read_file(mod.file, type = "MathProg")

class(mod)

## [1] "MP_data_from_file" "MILP"

res <- Rglpk_solve_LP(obj = mod$objective,

mat = mod$constraints[[1]],

dir = mod$constraints[[2]],

rhs = mod$constraints[[3]], bounds = mod

无约束的非线性规划

f(x1, x2) = (1 􀀀 x1)2 + 100(x2 􀀀 x21)2 .

Rosenbrock 香蕉函数

函数参数

par：表示各变量的初始值，通过向量传入

fn：目标函数

gr：梯度函数

method：优化方法，通过字符串的形式传入

control：一个列表，包含优化中的各种设置，包括最大迭代次

数maxit、绝对收敛容忍度abstol、相对收敛容忍度reltol 等

hessian：表示是否返回海塞矩阵，默认是FALSE

Nelder-Mead

一种单纯型法

非常稳健，效率也不低

CG

一种共轭梯度法

最速下降

容易陷入局部最优

BFGS

一种拟牛顿法，也称为变尺度法

具备牛顿法搜索快的特性的基础上又能有效地搜索全局最优

解

L-BFGS-B 是对该方法的一个优化，能够在优化的同时增加

箱型的约束条件

代码示例

obj.rosenbrock <- function(x) {

x1 <- x[1]

x2 <- x[2]

100 * (x2 - x1 * x1)^2 + (1 - x1)^2

}

optim(par = c(0, 3), fn = obj.rosenbrock)

## $par

## [1] 1.000300 1.000562

##

## $value

## [1] 2.301092e-07

##

## $counts

## function gradient

## 83 NA

##

## $convergence

## [1] 0

带约束的非线性规划

constrOptim 函数

f 表示目标函数

grad 表示梯度函数

theta 表示初始值

ui 和ci 表示线性约束条件

所有不等式的方向都是“>=”

第三方包

alabama

Rsolnp

代码示例

library(alabama)

hin.rosenbrock <- function(x) {

x1 <- x[1]

x2 <- x[2]

c(-x1^2 - x2^2 + 4, x1 / x2 -2)

}

constrOptim.nl(par=c(1,0.3), fn=obj.rosenbrock,

gr=gr.rosenbrock, hin=hin.rosenbrock)

## $par

## [1] 0.5174369 0.2587185

##

## $value

## [1] 0.2410077

##

## $counts

## function gradient

## 365 56

##

## $convergence

## [1] 0

遗传算法简介

一种随机算法，基于进化论的思想，模拟物种的世代，优胜劣汰

步骤

随机模拟初始物种作为初始值

目标函数代表个体的优秀程度

通过轮盘法决定繁殖，优秀的个体容易找到配偶

基因交换

基因突变机制

世代更替

代码实现

gene.reproduce <- function(ind) {

pop.val <- sapply(ind, gene.obj)

pop.prob <- 1 / (pop.val - min(pop.val) + 1)^0.2

pop.prob <- pop.prob/sum(pop.prob) #构造轮盘概率

pop.idx <- sample(seq_along(ind), 2,

prob = pop.prob) #抽样选出亲代

ind.gene.len <- sample(1:13, 1)

ind.gene.i <- sample(1:14, ind.gene.len)

new1 <- ind[[pop.idx[1]]]

new2 <- ind[[pop.idx[2]]]

new1[ind.gene.i]<-ind[[pop.idx[2]]][ind.gene.i]

new2[ind.gene.i]<-ind[[pop.idx[1]]][ind.gene.i]

ind.new <- new1 # 染色体交换

# 选择存活的子代与淘汰最弱的亲代：

if (gene.obj(new2)<gene.obj(new1)) ind.new <- new2

ind.old.i <- which.min(pop.prob)

if (gene.obj(ind.new)<gene.obj(ind[[ind.old.i]]))

ind[[ind.old.i]] <- ind.new

return(ind)

}