不同R优化包（pracma/nloptr）带约束优化结果差异排查

数据生成 #

首先生成带噪声的模拟数据，真实参数为A=1.34, B=0.5673, C=6.356, D=-1.234：

set.seed(1)
A = 1.34; B = 0.5673; C = 6.356; D = -1.234
x = seq(0.5, 20, length.out = 500)
y = A + B * x + C * x^2 + D * log(x) + runif(500, 0, 3)
Rand_Data = data.frame(x = x, Val = y)

优化问题定义 #

目标是估计参数a,b,c,d，最小化观测值与拟合值的绝对差之和（L1损失），同时满足两个约束：

• 约束1：第100个拟合值小于120，即a + b*x[100] + c*x[100]^2 + d*log(x[100]) < 120
• 约束2：拟合值总和与观测值总和的差值绝对值小于0.001，即|sum(拟合值) - sum(Rand_Data$Val)| < 0.001

拟合函数定义：

MyFN = function(a, b, c, d) return(a + b * Rand_Data$x + c * Rand_Data$x^2 + d * log(Rand_Data$x))

现有优化结果 #

使用pracma包的实现 #

通过fmincon正确定义线性约束和非线性不等式约束：

library(pracma)
eps <- 1e-4
X <- cbind(1, x, x^2, log(x))
f <- function(theta) {
  sum(abs(X %*% theta - y))
}
# 线性约束：X[100,]θ ≤ 120 - eps（对应约束1）
A <- X[100, , drop = FALSE]
b <- 120 - eps
# 非线性不等式约束：abs(sum(Xθ)-sum(y)) -1e-3 + eps ≤0（对应约束2）
hin <- function(theta) {
  abs(sum(X %*% theta) - sum(y)) - 1e-3 + eps
}
sol_pracma = fmincon(rep(0, 4), f, A = A, b = b, hin = hin)$par
sol_pracma # 输出：774.972960   318.128595     4.354805 -1798.378655

使用nloptr包的实现 #

原代码错误地将约束1设为等式约束，导致结果偏离：

library(nloptr)
Hx <- function(theta) {
    X[100, , drop = FALSE] %*% theta - (120 - eps)
}
print(nloptr(rep(0, 4), f, eval_g_ineq = hin, eval_g_eq = Hx, opts = list("algorithm"="NLOPT_LN_COBYLA", "xtol_rel"=1.0e-8))) 

# 输出参数：0.3596822 -0.8529221 6.461235 0.03229311

差异原因分析 #

1. 约束定义错误：nloptr代码中把约束1（不等式约束MyFN[100]<120）误设为等式约束（eval_g_eq），强行要求第100个拟合值等于120-eps，完全违背原约束要求；而pracma中正确将其定义为线性不等式约束。
2. 算法特性差异：pracma的fmincon默认使用SQP（序列二次规划）算法，适合处理带线性/非线性约束的优化问题；nloptr选用的NLOPT_LN_COBYLA是导数自由的非线性约束优化算法，对初始值和约束定义敏感度更高，加上约束错误，导致结果严重偏离。
3. 约束处理逻辑不一致：pracma中hin函数要求返回值≤0，对应约束2的差值绝对值小于0.001；nloptr原代码虽复用了hin函数，但约束1的错误定义直接改变了优化问题的可行域。

解决方法 #

1. 修正nloptr的约束定义 #

将两个约束都正确定义为不等式约束，确保可行域符合原问题要求：

library(nloptr)
eps <- 1e-4
X <- cbind(1, x, x^2, log(x))
y <- Rand_Data$Val

# 目标函数
f <- function(theta) {
  sum(abs(X %*% theta - y))
}

# 定义所有不等式约束：返回值需 ≤0
g_ineq <- function(theta) {
  c(
    # 约束1：第100个拟合值 <120 → X[100,]θ - (120 - eps) ≤0
    X[100, , drop = FALSE] %*% theta - (120 - eps),
    # 约束2：总和差值绝对值 <0.001 → abs(sum(Xθ)-sum(y)) - (1e-3 - eps) ≤0
    abs(sum(X %*% theta) - sum(y)) - (1e-3 - eps)
  )
}

# 运行优化，调整迭代次数确保收敛
sol_nloptr <- nloptr(
  x0 = rep(0, 4),
  eval_f = f,
  eval_g_ineq = g_ineq,
  opts = list(
    "algorithm" = "NLOPT_LN_COBYLA",
    "xtol_rel" = 1.0e-8,
    "maxeval" = 10000
  )
)

print(sol_nloptr$solution)

2. 验证约束满足情况 #

优化后需检查两个约束是否都满足：

# 检查pracma结果
theta_pracma <- sol_pracma
cat("pracma结果-第100个拟合值：", X[100,] %*% theta_pracma, "\n")
cat("pracma结果-总和差值绝对值：", abs(sum(X %*% theta_pracma) - sum(y)), "\n")

# 检查修正后nloptr结果
theta_nloptr <- sol_nloptr$solution
cat("nloptr结果-第100个拟合值：", X[100,] %*% theta_nloptr, "\n")
cat("nloptr结果-总和差值绝对值：", abs(sum(X %*% theta_nloptr) - sum(y)), "\n")

3. 选择适配的算法 #

对于L1损失的约束优化，可根据情况选择：

• 若问题存在线性约束，优先使用pracma的fmincon（SQP算法），收敛性更稳定；
• 使用nloptr时，可尝试NLOPT_LN_SBPLX算法，或调整COBYLA的maxeval参数确保收敛到可行域内的最优解。

内容的提问来源于stack exchange，提问作者Daniel Lobo