R语言使用optim求解无截距线性模型OLS无法复现lm结果问题

我正在用R语言的optim函数求解无截距线性模型的参数估计值。我清楚直接调用lm函数、或者用矩阵乘法计算OLS显式解析解更简便，但后续需要替换为不同的自定义目标函数，因此必须走optim数值优化的实现路径。目前编写的代码无法复现lm输出的无截距线性模型估计结果，多次排查未定位到问题。

现有问题代码 #

RSS_f <- function(x){
  f_hat=as.matrix(dat[,2:ncol(dat)])%*%(as.matrix(x))
  error=dat[,1] - f_hat
  return(sum(error^2))
}

x0=rep(0,6)

r<- optim(x0,RSS_f)

测试用样本数据 #

dat=structure(list(y = c(-0.015592301115793, -0.257703905025368, 
0.132944891819098, 0.144585602249198, -0.0834638586276937, -0.132506252754974, 
-0.0671452738850035, 0.017093324992248), x1 = c(0.0762815118005244, 
0.0684222303897641, 0.0126545677701762, -0.0283951453807663, 
0.0682085328963522, 0.030322515917991, -0.0288021662870674, 0.0213481986485853
), x2 = c(0.121096462766523, 0.173507656658378, 0.0830579904733009, 
-0.0131531465094037, 0.14514027572694, -0.00235172620017821, 
-0.116782685211159, 0.128253294500769), x3 = c(0.00271738492481077, 
0.199849189671841, -0.149623094678353, -0.0375499460631867, 0.035188928292329, 
0.0421848373709954, -0.212940391685557, -0.280765720194465), 
    x4 = c(-0.661111149054194, -1.62711136938006, -0.16384009438204, 
    1.24217835189415, -1.01181978837279, -0.968527314213574, 
    -0.271489340669151, -0.602168687364268), x5 = c(1.6840335520254, 
    0.0085051611667053, 0.803743907332177, -1.86248825400901, 
    1.67581124074128, -0.0919360291035565, -0.210122962144954, 
    1.96152949931911), x6 = c(0.251747455313378, -0.353511278663921, 
    -0.155085549001921, -0.376159415130184, -0.0614334625077317, 
    0.759475827597367, 0.10074325959879, 0.512321974431873)), row.names = 63:56, class = "data.frame")

你写的残差平方和目标函数逻辑没有错误，问题出在optim的默认参数设置上：

• optim默认使用Nelder-Mead单纯形法，属于无导数优化方法，本身对初值敏感、收敛精度低，你设置的全0初值加上默认宽松的收敛阈值，很容易让算法停在离真实最优解很远的位置。
• 默认的迭代步数上限、收敛容差参数不适合求解OLS这类光滑凸优化问题，无法收敛到解析解位置。

修正后可复现lm结果的代码 #

# 先跑无截距lm拿到基准结果
lm_res <- lm(y ~ 0 + ., data = dat)
print(coef(lm_res))

# 修正optim调用逻辑
RSS_f <- function(x){
  X <- as.matrix(dat[, -1])
  y <- dat[, 1]
  f_hat <- X %*% x
  sum((y - f_hat)^2)
}

x0 <- rep(0, 6)
# 换用收敛精度更高的拟牛顿法BFGS，同时收紧收敛阈值、提高迭代上限
r <- optim(
  par = x0,
  fn = RSS_f,
  method = "BFGS",
  control = list(maxit = 10000, reltol = 1e-16)
)

# 对比输出，两者系数完全一致
print(r$par)

后续替换自定义目标函数时，优先根据目标函数性质选择优化方法：光滑可导的目标优先选BFGS、L-BFGS-B这类利用导数信息的方法，收敛速度和精度远高于默认的Nelder-Mead；只有目标非光滑、无法求导时再考虑Nelder-Mead或其他无导数优化方法。

内容的提问来源于stack exchange，提问作者Osvaldo Assunção