使用 nloptr 与 Excel 在 R 中进行非线性优化
Non Linear Optimization in R with nloptr vs Excel
我似乎找不到我的特定问题的答案。我正在尝试最小化实值向量和模型线性组合之间的平均绝对误差,如下所示:
library(nloptr)
df <- data.frame(
real = c(24.2418, 21.7374, 7.203, 115.0233, 16.632, 5.4644, 27.8917, 0.5904, 0.633, 105.3656, 110.0564, 122.9399, 23.0418, 4.2186, 109.5453),
m_1 = c(17.7790979854662, 32.93639375, 6.94294375000066, 98.909065625, 80.1068562499999, 11.656556250002, 39.868921875, 0.859157480988586, 0.612625482376216, 112.580383416359, 151.896765176957, 155.81521460987, 7.3257, 4.1268, 41.5711703205879),
m_2 = c(25.5062900474607, 32.7709877317137, 6.94317649489279, 98.898593448344, 80.1114267521597, 11.6563450007055, 39.8691722409057, 0.907260912997819, 0.602795676399197, 114.183526809465, 139.51052151724, 149.993624420536, 6.85002142907728, 4.66668862149305, 70.7527906311631),
m_3 = c(27.1495912199999, 40.2339353399999, 7.10339542, 87.1444967133334, 58.4563384933332, 11.1378198366666, 37.6030141333333, 0.852288459999999, 0.681724560000002, 100.101136186666, 118.536525109999, 136.923264319999, 5.64763034333333, 4.8659515, 70.12675835),
m_4 = c(25.511590625, 32.9363937499999, 7.00050769504031, 98.3660974929738, 80.10685625, 11.65655625, 39.868921875, 0.665580984791989, 0.498756215272925, 85.791042265746, 135.619508469251, 140.946144588455, 5.05824305930683, 3.25333636845094, 22.0908752674237),
m_5 = c(25.6118152340655, 34.5646719234769, 6.82416840383483, 91.5582383465651, 84.4603847826215, 11.3405620277701, 40.7906062254877, 0.908706704665592, 0.602817399156822, 114.326905157898, 139.595783699511, 150.046375909198, 6.8511793011574, 4.6622942290559, 56.2753212961812),
m_6 = c(21.9868574376585, 44.3147731773466, 6.38315486686481, 100.303757097094, 9.13921010739697, 7.83817900918309, 31.5458855316741, 1.09960505333834, 0.817751834425678, 101.110814846224, 145.55847538105, 142.82362305075, 7.61732986965459, 4.6774198307473, 67.5821464371521)
)
best_dist <- function(x) {
output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] +
df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
mean(abs(output - df$real))
}
restriction <- function(x) sum(x) - 1
nloptr(
x0 = rep(1 / 6, 6),
eval_f = best_dist,
lb = rep(0, 6),
ub = rep(1, 6),
eval_g_eq = restriction,
opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)
如您所见,我正在使用 nloptr
包。上面的代码为 objective 函数产生了 14.85 的非最优结果,参数都是初始参数。您可能将初始参数更改为其他一些向量,但仍然无法获得最优解。
但是,使用 excel 解算器可以轻松得到 objective 函数的 10.77 和 (0, 0, .15, 0, 0, .85) 参数。
我试过使用带有渐变的算法,但我似乎无法获得正确的语法。这是我的另一个尝试。
gradient <- function(x) {
output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] +
df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
err <- output - df$real
c(
- sum(sign(err) * df$m_1),
- sum(sign(err) * df$m_2),
- sum(sign(err) * df$m_3),
- sum(sign(err) * df$m_4),
- sum(sign(err) * df$m_5),
- sum(sign(err) * df$m_6)
)
}
nloptr(
x0 = runif(6),
eval_f = best_dist,
eval_grad_f = gradient,
lb = rep(0, 6),
ub = rep(1, 6),
eval_g_eq = restriction,
opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)
这类似于 LAD 回归,带有边约束。显示了 LAD 回归的 LP(线性规划)公式 here。使用您对数据的命名 (real[i], m[i,j]),我们有:
min sum(i, | real[i] - sum(j,m[i,j]*x[j]) | )
subject to
sum(j, x[j]) = 1
0 <= x[j] <= 1
这可以线性化(使用变量拆分)如下:
min sum(i, r1[i]+r2[i])
subject to
r1[i]-r2[i] = real[i] - sum(j,m[i,j]*x[j])
sum(j, x[j]) = 1
bounds:
r1[i], r2[i] >= 0 (positive variables)
0 <= x[j] <= 1
这是一个纯 LP,可以用任何 LP 求解器求解。
我似乎找不到我的特定问题的答案。我正在尝试最小化实值向量和模型线性组合之间的平均绝对误差,如下所示:
library(nloptr)
df <- data.frame(
real = c(24.2418, 21.7374, 7.203, 115.0233, 16.632, 5.4644, 27.8917, 0.5904, 0.633, 105.3656, 110.0564, 122.9399, 23.0418, 4.2186, 109.5453),
m_1 = c(17.7790979854662, 32.93639375, 6.94294375000066, 98.909065625, 80.1068562499999, 11.656556250002, 39.868921875, 0.859157480988586, 0.612625482376216, 112.580383416359, 151.896765176957, 155.81521460987, 7.3257, 4.1268, 41.5711703205879),
m_2 = c(25.5062900474607, 32.7709877317137, 6.94317649489279, 98.898593448344, 80.1114267521597, 11.6563450007055, 39.8691722409057, 0.907260912997819, 0.602795676399197, 114.183526809465, 139.51052151724, 149.993624420536, 6.85002142907728, 4.66668862149305, 70.7527906311631),
m_3 = c(27.1495912199999, 40.2339353399999, 7.10339542, 87.1444967133334, 58.4563384933332, 11.1378198366666, 37.6030141333333, 0.852288459999999, 0.681724560000002, 100.101136186666, 118.536525109999, 136.923264319999, 5.64763034333333, 4.8659515, 70.12675835),
m_4 = c(25.511590625, 32.9363937499999, 7.00050769504031, 98.3660974929738, 80.10685625, 11.65655625, 39.868921875, 0.665580984791989, 0.498756215272925, 85.791042265746, 135.619508469251, 140.946144588455, 5.05824305930683, 3.25333636845094, 22.0908752674237),
m_5 = c(25.6118152340655, 34.5646719234769, 6.82416840383483, 91.5582383465651, 84.4603847826215, 11.3405620277701, 40.7906062254877, 0.908706704665592, 0.602817399156822, 114.326905157898, 139.595783699511, 150.046375909198, 6.8511793011574, 4.6622942290559, 56.2753212961812),
m_6 = c(21.9868574376585, 44.3147731773466, 6.38315486686481, 100.303757097094, 9.13921010739697, 7.83817900918309, 31.5458855316741, 1.09960505333834, 0.817751834425678, 101.110814846224, 145.55847538105, 142.82362305075, 7.61732986965459, 4.6774198307473, 67.5821464371521)
)
best_dist <- function(x) {
output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] +
df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
mean(abs(output - df$real))
}
restriction <- function(x) sum(x) - 1
nloptr(
x0 = rep(1 / 6, 6),
eval_f = best_dist,
lb = rep(0, 6),
ub = rep(1, 6),
eval_g_eq = restriction,
opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)
如您所见,我正在使用 nloptr
包。上面的代码为 objective 函数产生了 14.85 的非最优结果,参数都是初始参数。您可能将初始参数更改为其他一些向量,但仍然无法获得最优解。
但是,使用 excel 解算器可以轻松得到 objective 函数的 10.77 和 (0, 0, .15, 0, 0, .85) 参数。
我试过使用带有渐变的算法,但我似乎无法获得正确的语法。这是我的另一个尝试。
gradient <- function(x) {
output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] +
df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
err <- output - df$real
c(
- sum(sign(err) * df$m_1),
- sum(sign(err) * df$m_2),
- sum(sign(err) * df$m_3),
- sum(sign(err) * df$m_4),
- sum(sign(err) * df$m_5),
- sum(sign(err) * df$m_6)
)
}
nloptr(
x0 = runif(6),
eval_f = best_dist,
eval_grad_f = gradient,
lb = rep(0, 6),
ub = rep(1, 6),
eval_g_eq = restriction,
opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)
这类似于 LAD 回归,带有边约束。显示了 LAD 回归的 LP(线性规划)公式 here。使用您对数据的命名 (real[i], m[i,j]),我们有:
min sum(i, | real[i] - sum(j,m[i,j]*x[j]) | )
subject to
sum(j, x[j]) = 1
0 <= x[j] <= 1
这可以线性化(使用变量拆分)如下:
min sum(i, r1[i]+r2[i])
subject to
r1[i]-r2[i] = real[i] - sum(j,m[i,j]*x[j])
sum(j, x[j]) = 1
bounds:
r1[i], r2[i] >= 0 (positive variables)
0 <= x[j] <= 1
这是一个纯 LP,可以用任何 LP 求解器求解。