R优化与Scipy优化之间的差异：Nelder-Mead

Question

我写了一个脚本，我认为它应该在Python和R中产生相同的结果，但它们产生了非常不同的答案。每个都试图使用Nelder-Mead最小化偏差来拟合模拟数据。总的来说，R中的优化表现要好得多。我做错了什么吗？R和SciPy中实现的算法不同吗？

Python的结果：

>>> res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')

 final_simplex: (array([[-0.21483287, -1.        , -0.4645897 , -4.65108495],
       [-0.21483909, -1.        , -0.4645915 , -4.65114839],
       [-0.21485426, -1.        , -0.46457789, -4.65107337],
       [-0.21483727, -1.        , -0.46459331, -4.65115965],
       [-0.21484398, -1.        , -0.46457725, -4.65099805]]), array([107.46037865, 107.46037868, 107.4603787 , 107.46037875,
       107.46037875]))
           fun: 107.4603786452194
       message: 'Optimization terminated successfully.'
          nfev: 349
           nit: 197
        status: 0
       success: True
             x: array([-0.21483287, -1.        , -0.4645897 , -4.65108495])

R结果：

> res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

$par
[1] 0.2641022 1.0000000 0.2086496 3.6688737

$value
[1] 110.4249

$counts
function gradient 
     329       NA 

$convergence
[1] 0

$message
NULL

我已检查过我的代码，据我所知，这似乎是由于优化和最小化之间的某些差异导致的，因为我试图最小化的函数（即choiceProbDev）在每个函数中的操作方式都相同（除了输出外，我还检查了函数中每个步骤的等价性）。例如：

Python的choiceProbDev函数：

>>> choiceProbDev(np.array([0.5, 0.5, 0.5, 3]), stim, dflt, dat, N)
143.31438613033876

R选择概率开发：

> choiceProbDev(c(0.5, 0.5, 0.5, 3), stim, dflt, dat, N)
[1] 143.3144

我也尝试了调整每个优化函数的公差水平，但我不确定这两个参数如何匹配。无论如何，目前我的摆弄并没有使它们达成一致。以下是每个代码的完整内容。

Python:

# load modules
import math
import numpy as np
from scipy.optimize import minimize
from scipy.stats import binom

# initialize values
dflt = 0.5
N = 1

# set the known parameter values for generating data
b = 0.1
w1 = 0.75
w2 = 0.25
t = 7

theta = [b, w1, w2, t]

# generate stimuli
stim = np.array(np.meshgrid(np.arange(0, 1.1, 0.1),
                            np.arange(0, 1.1, 0.1))).T.reshape(-1,2)

# starting values
sparams = [-0.5, -0.5, -0.5, 4]


# generate probability of accepting proposal
def choiceProb(stim, dflt, theta):

    utilProp = theta[0] + theta[1]*stim[:,0] + theta[2]*stim[:,1]  # proposal utility
    utilDflt = theta[1]*dflt + theta[2]*dflt  # default utility
    choiceProb = 1/(1 + np.exp(-1*theta[3]*(utilProp - utilDflt)))  # probability of choosing proposal

    return choiceProb

# calculate deviance
def choiceProbDev(theta, stim, dflt, dat, N):

    # restrict b, w1, w2 weights to between -1 and 1
    if any([x > 1 or x < -1 for x in theta[:-1]]):
        return 10000

    # initialize
    nDat = dat.shape[0]
    dev = np.array([np.nan]*nDat)

    # for each trial, calculate deviance
    p = choiceProb(stim, dflt, theta)
    lk = binom.pmf(dat, N, p)

    for i in range(nDat):
        if math.isclose(lk[i], 0):
            dev[i] = 10000
        else:
            dev[i] = -2*np.log(lk[i])

    return np.sum(dev)


# simulate data
probs = choiceProb(stim, dflt, theta)

# randomly generated data based on the calculated probabilities
# dat = np.random.binomial(1, probs, probs.shape[0])
dat = np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
       0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
       0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
       0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
       0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

# fit model
res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')

R:

library(tidyverse)

# initialize values
dflt <- 0.5
N <- 1

# set the known parameter values for generating data
b <- 0.1
w1 <- 0.75
w2 <- 0.25
t <- 7

theta <- c(b, w1, w2, t)

# generate stimuli
stim <- expand.grid(seq(0, 1, 0.1),
                    seq(0, 1, 0.1)) %>%
  dplyr::arrange(Var1, Var2)

# starting values
sparams <- c(-0.5, -0.5, -0.5, 4)

# generate probability of accepting proposal
choiceProb <- function(stim, dflt, theta){
  utilProp <- theta[1] + theta[2]*stim[,1] + theta[3]*stim[,2]  # proposal utility
  utilDflt <- theta[2]*dflt + theta[3]*dflt  # default utility
  choiceProb <- 1/(1 + exp(-1*theta[4]*(utilProp - utilDflt)))  # probability of choosing proposal
  return(choiceProb)
}

# calculate deviance
choiceProbDev <- function(theta, stim, dflt, dat, N){
  # restrict b, w1, w2 weights to between -1 and 1
  if (any(theta[1:3] > 1 | theta[1:3] < -1)){
    return(10000)
  }

  # initialize
  nDat <- length(dat)
  dev <- rep(NA, nDat)

  # for each trial, calculate deviance
  p <- choiceProb(stim, dflt, theta)
  lk <- dbinom(dat, N, p)

  for (i in 1:nDat){
    if (dplyr::near(lk[i], 0)){
      dev[i] <- 10000
    } else {
      dev[i] <- -2*log(lk[i])
    }
  }
  return(sum(dev))
}

# simulate data
probs <- choiceProb(stim, dflt, theta)

# same data as in python script
dat <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
         0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
         0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

# fit model
res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

更新:

在每个迭代中打印估计值后，我现在认为差异可能源于每个算法所采取的“步长”不同。Scipy似乎比optim采取更小的步长(并且初始方向也不同)。我还没有想出如何调整这一点。

Python:

>>> res = minimize(choiceProbDev, sparams, (stim, dflt, dat, N), method='Nelder-Mead')
[-0.5 -0.5 -0.5  4. ]
[-0.525 -0.5   -0.5    4.   ]
[-0.5   -0.525 -0.5    4.   ]
[-0.5   -0.5   -0.525  4.   ]
[-0.5 -0.5 -0.5  4.2]
[-0.5125 -0.5125 -0.5125  3.8   ]
...

R:

> res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N, method="Nelder-Mead")
[1] -0.5 -0.5 -0.5  4.0
[1] -0.1 -0.5 -0.5  4.0
[1] -0.5 -0.1 -0.5  4.0
[1] -0.5 -0.5 -0.1  4.0
[1] -0.5 -0.5 -0.5  4.4
[1] -0.3 -0.3 -0.3  3.6
...

- YTD

---

这是一个非常复杂的使用案例，不是简单的“optim”。我建议打印出数据处理过程的每个步骤，并检查两个脚本中的每个数据部分是否匹配。这些步骤中的任何一个都可能是问题所在。其中之一检查：“binom.pmf”和“dbinom”。 - Parfait

@Parfait -- 谢谢！我编辑了我的问题，希望能够突出和准确定位问题。我已经编辑了我的代码，更清晰地指定了正在使用的库。我逐行检查了两个代码，一切看起来都是等价的（例如，我正在最小化的choiceProbDev函数在两个实现中产生相同的结果）。 - YTD

尝试在R的optim control list中玩弄args，以与Python的minimize(method='Nelder-Mead')默认值匹配。我认为这是由于默认值不同造成的。 - Parfait

optim 中的 Nelder-Mead 方法并不是 R 中可用的最佳或最准确的实现。您可以尝试使用 dfoptim 包中的 nmk[b] 或 adagio 中的 neldermead[b]，或者在 pracma 中尝试自适应版本 anms 等其他方法。这些实现并没有太大的区别，但在精度和效率方面可能会有显著差异，特别是如果存在多个最小值时。 - Hans W.

Answer 1

这并不完全是关于“优化器差异”的答案，但我想在这里贡献一些对优化问题的探索。以下是几个要点：

• 表面很平滑，所以基于导数的优化器可能效果更好（即使没有显式编写梯度函数，例如退回到有限差分逼近 - 如果有梯度函数则会更好）
• 该表面是对称的，因此具有多个最优解（明显有两个），但它不是高度多峰或粗糙，因此我认为随机全局优化器不值得麻烦
• 对于不太高维或计算成本较高的优化问题，可视化全局表面以了解情况是可行的。
• 对于带约束的优化，通常最好 使用显式处理约束的优化器，或者 将参数比例更改为无约束比例

这是整个表面的图片：

红色轮廓是对数似然等于（110、115、120）的轮廓（我能得到的最佳拟合是LL = 105.7）。 最佳点位于第二列第三行（通过L-BFGS-B实现）和第五列第四行（真实参数值）。 （我还没有检查目标函数以查看对称性来自哪里，但我认为这可能很清楚。）Python的Nelder-Mead和R的Nelder-Mead表现几乎相同。

---

参数和问题设置

## initialize values
dflt <- 0.5; N <- 1
# set the known parameter values for generating data
b <- 0.1; w1 <- 0.75; w2 <- 0.25; t <- 7
theta <- c(b, w1, w2, t)
# generate stimuli
stim <- expand.grid(seq(0, 1, 0.1), seq(0, 1, 0.1))
# starting values
sparams <- c(-0.5, -0.5, -0.5, 4)
# same data as in python script
dat <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,
         0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
         0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
         0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)

目标函数

尽可能使用内置函数（例如plogis()，dbinom(...,log=TRUE)）。

# generate probability of accepting proposal
choiceProb <- function(stim, dflt, theta){
    utilProp <- theta[1] + theta[2]*stim[,1] + theta[3]*stim[,2]  # proposal utility
    utilDflt <- theta[2]*dflt + theta[3]*dflt  # default utility
    choiceProb <- plogis(theta[4]*(utilProp - utilDflt))  # probability of choosing proposal
    return(choiceProb)
}
# calculate deviance
choiceProbDev <- function(theta, stim, dflt, dat, N){
  # restrict b, w1, w2 weights to between -1 and 1
    if (any(theta[1:3] > 1 | theta[1:3] < -1)){
        return(10000)
    }
    ## for each trial, calculate deviance
    p <-  choiceProb(stim, dflt, theta)
    lk <-  dbinom(dat, N, p, log=TRUE)
    return(sum(-2*lk))
}
# simulate data
probs <- choiceProb(stim, dflt, theta)

模型拟合

# fit model
res <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")
## try derivative-based, box-constrained optimizer
res3 <- optim(sparams, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
              lower=c(-1,-1,-1,-Inf), upper=c(1,1,1,Inf),
             method="L-BFGS-B")

py_coefs <- c(-0.21483287,  -0.4645897 , -1, -4.65108495) ## transposed?
true_coefs <- c(0.1, 0.25, 0.75, 7)  ## transposed?
## start from python coeffs
res2 <- optim(py_coefs, choiceProbDev, stim=stim, dflt=dflt, dat=dat, N=N,
             method="Nelder-Mead")

探索对数似然表面

cc <- expand.grid(seq(-1,1,length.out=51),
                  seq(-1,1,length.out=6),
                  seq(-1,1,length.out=6),
                  seq(-8,8,length.out=51))
## utility function for combining parameter values
bfun <- function(x,grid_vars=c("Var2","Var3"),grid_rng=seq(-1,1,length.out=6),
                 type=NULL) {
    if (is.list(x)) {
        v <- c(x$par,x$value)
    } else if (length(x)==4) {
        v <- c(x,NA)
    }
    res <- as.data.frame(rbind(setNames(v,c(paste0("Var",1:4),"z"))))
    for (v in grid_vars)
        res[,v] <- grid_rng[which.min(abs(grid_rng-res[,v]))]
    if (!is.null(type)) res$type <- type
    res
}

resdat <- rbind(bfun(res3,type="R_LBFGSB"),
                bfun(res,type="R_NM"),
                bfun(py_coefs,type="Py_NM"),
                bfun(true_coefs,type="true"))

cc$z <- apply(cc,1,function(x) choiceProbDev(unlist(x), dat=dat, stim=stim, dflt=dflt, N=N))
library(ggplot2)
library(viridisLite)
ggplot(cc,aes(Var1,Var4,fill=z))+
    geom_tile()+
    facet_grid(Var2~Var3,labeller=label_both)+
    scale_fill_viridis_c()+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0))+
    theme(panel.spacing=grid::unit(0,"lines"))+
    geom_contour(aes(z=z),colour="red",breaks=seq(105,120,by=5),alpha=0.5)+
    geom_point(data=resdat,aes(colour=type,shape=type))+
    scale_colour_brewer(palette="Set1")

ggsave("liksurf.png",width=8,height=8)

- Ben Bolker

Answer 2

"'Nelder-Mead'一直是一个问题优化方法，并且在optim中的编码已经过时。我们将尝试R包中提供的其他三种实现。为了避免其他参数，让我们将函数fn定义为："

fn <- function(theta)
        choiceProbDev(theta, stim=stim, dflt=dflt, dat=dat, N=N)

然后，求解器dfoptim::nmk()、adagio::neldermead()和pracma::anms()将都返回相同的最小值xmin = 105.7843，但在不同的位置，例如

dfoptim::nmk(sparams, fn)
## $par
## [1] 0.1274937 0.6671353 0.1919542 8.1731618
## $value
## [1] 105.7843

这些是真正的局部极小值，而例如 Python 解决方案 107.46038 在 c(-0.21483287,-1.0,-0.4645897,-4.65108495) 上不是。显然，您的问题数据不足以拟合模型。

您可以尝试全局优化器，在一定范围内可能会找到全局最优解。对我来说，所有局部最小值具有相同的最小值。

- Hans W.

---

感谢您提供全局优化器的建议。同时，也感谢您提供了其他可能更好的R包的有益建议。如果optim已经如此过时，那么它在这种情况下表现得比scipy.optimize好得多是很有趣的。难道scipy还更过时吗？我想我的问题仍然存在：为什么scipy的结果与所有列出的R优化器相比如此不同且差距如此之大？我想了解差异的原因，以及是否有办法提高其性能/使其与R更接近。 - YTD

此外，是否有另一个Python模块可以更好地执行Nelder-Mead算法？每当我在搜索Python中的优化函数时，似乎绝大部分结果都指向scipy.optimize。毫无疑问，Python一定有一种可比较的算法，至少可以接近R的算法。 - YTD

为什么您认为optim的性能比scipy.minimize要好得多？ optim返回110.4249作为最小值，而minimize默认选项下返回107.4604。两者都不是真正的局部极小值。 SciPy源代码只提到Nelder和Mead的原始文章以及1996年的概述文章。将自适应的Nelder-Mead过程（大约20 LoC）从R（或其Matlab版本）转换为Python并不太困难。 - Hans W.

请查看PyPi项目nelder-mead，最新版本为2018年10月，这可能是为Python社区提供新实现的尝试。使用Python3可以通过pip install nelder-mead进行安装。（免责声明：我没有尝试过。） - Hans W.

啊，你说得对……我一直关注于optim始终返回比原始输入更接近的参数值这一事实，以至于忽略了从scipy.optimize获得的最小值实际上更低的事实。糟糕。谢谢你提供的PyPi提示。 - YTD