在这种情况下,为什么 L_BFGS_B 优化会跳到可行解决方案的极端范围?

Why does L_BFGS_B optimization skip to extreme range of viable solutions in this instance?

我很感激这是一个非常具体的问题!

为了帮助解释:我正在探索使用线性优化器来演示函数曲面中的尖锐 'cliffs' 会导致非最佳解决方案。 R 中的可重现代码如下:

library(glmnet)
library(mice)


# Load data
df <- read.csv(paste0('https://raw.githubusercontent.com/jbrownlee/Datasets',
                      '/master/pima-indians-diabetes.data.csv'), header = F)

colnames(df) <- c('Pregnancies', 'Glucose', 'BloodPressure', 'SkinThickness',
                  'Insulin', 'BMI', 'DiabetesPedigreeFunction', 'Age', 'Outcome')


set.seed(40)

# Impute 0 (missing) values for columns 2 through 8 (Glucose - Age)
df[2:8] <- lapply(df[2:8], function(x) replace(x, x %in% 0, NA))
micedf <- mice(df)
df <- complete(micedf)

# Create train/test split
sample_size <- floor(0.75 * nrow(df))
train_index <- sample(seq_len(nrow(df)), size = sample_size)
train <- df[train_index,]
test  <- df[-train_index,]

# Generate model matrix format for glmnet
x <- as.matrix(train[,1:8])
y <- train$Outcome

# Fitting function
GLM_tune <- function(alpha) {
    set.seed(40)
    cvglmnet <- glmnet::cv.glmnet(x, y, nfolds = 5, family = "binomial",
                                  alpha = alpha, type.measure = "auc",
                                  parallel = F)

    return (cvglmnet$cvm[cvglmnet$lambda == cvglmnet$lambda.1se])    }

现在,如果我输入一个介于 0 和 1 之间的值,如下所示:

optim(par = 0.9, fn = GLM_tune, lower = 0, upper = 1, 
      control = list(fnscale = -1, trace=3), method = c("L-BFGS-B"))

# >> $par = 0.86

优化器爬升到局部最大值 - 我已经通过使用以下方法探索整个表面区域进行了测试:

surf <- data.frame(alpha = 0, auc = 0)   
for (a in seq(from=0, to=1000)) {
    surf[a+1,1] <- a/1000
    surf[a+1,2] <- GLM_tune(a/1000)
}

library(ggplot2)
ggplot() +
    geom_point(data=surf, size = 1.2, color = "black", aes(alpha, auc))

但是,当我将 alpha = 1 设置为起点时,算法在第二次迭代时转换为 alpha = 0,然后作为 'final' 解决方案退出:

optim(par = 1, fn = GLM_tune, lower = 0, upper = 1, 
      control = list(fnscale = -1, trace=3), method = c("L-BFGS-B"))

# >> $par = 0

为什么会这样?显然我不完全理解该算法,但我假设 optim 函数中的默认步进为 0.001(请参阅 ndeps)——那么为什么下一步会走向相反的极端?我是否缺少应针对这些问题设置的关键参数?

从 objective 函数的路径可以清楚地看出它有很多局部最大值,因此,像 "L-BFGS-B" 这样的基于梯度的优化算法 不是 适合寻找全局最大值。

此外,我的R(3.6),

optim(par = 1, fn = GLM_tune, lower = 0, upper = 1, 
      control = list(fnscale = -1, trace=3), method = c("L-BFGS-B"))$par
## [1] 1

returns 1 而不是您指出的 0

要理解为什么收敛到1,我们可以看一下"L-BFGS-B"算法的优化路径。为此,我更喜欢使用 R 包 optimParallel https://CRAN.R-project.org/package=optimParallel。我是包的作者:

library("optimParallel")
cl <- makeCluster(2); setDefaultCluster(cl=cl)
clusterExport(cl, c("x", "y")) # export implicitly used values
optimParallel(par = 1, fn = GLM_tune,
              lower = 0, upper = 1, 
              control = list(fnscale = -1),
              parallel = list(optimParallel.loginfo=TRUE))$loginfo 
##      step       par1         fn        gr1 
## 1.0000000  1.0000000 -0.8215854  0.0000000 

我们看到 1 处的梯度是 0。因此,算法停在 1.

也就不足为奇了

我们可以用

检查近似梯度的计算
ndeps <- 0.001  # the default value
(GLM_tune(1) - GLM_tune(1-ndeps))/ndeps
## [1] 0

请注意,如果 1 不是上限,则 optim() 使用中心差分梯度近似。像

(GLM_tune(1+ndeps) - GLM_tune(1-ndeps))/(2*ndeps)