解决R中nls()拟合奇异梯度错误的初始值设置问询

尝试为520个用户的数据拟合指数模型 Value ~ A * exp(-k * x)，调用nls()时反复出现奇异梯度错误。推测该错误与初始值设置相关，但无法找到适配所有用户的通用初始值，需要获取全部用户的拟合系数A和k用于后续分析。

数据示例 #

dput(head(Mfrq.df.2))
structure(list(User.ID = c("37593", "38643", "49433", "60403", 
"70923", "85363"), V1 = c(9L, 3L, 4L, 80L, 19L, 0L), V2 = c(10L, 
0L, 29L, 113L, 21L, 1L), V3 = c(5L, 2L, 17L, 77L, 7L, 2L), V4 = c(2L, 
2L, 16L, 47L, 4L, 3L), V5 = c(2L, 10L, 16L, 40L, 1L, 8L), V6 = c(4L, 
0L, 9L, 22L, 1L, 7L), V7 = c(6L, 8L, 9L, 8L, 0L, 6L), V8 = c(2L, 
17L, 16L, 24L, 2L, 1L), V9 = c(3L, 20L, 7L, 30L, 0L, 4L), V10 = c(2L, 
11L, 5L, 11L, 2L, 3L)), row.names = c(NA, 6L), class = "data.frame")

已尝试的两种实现代码 #

方法一：分组拟合 #

# Way I
x=1:10
Mfrq.df.2_long <- pivot_longer(Mfrq.df.2, matches("V\\d{1,2}"), names_to = NULL, values_to = "Value")

Mfrq.df.2_long %>% 
  group_by(User.ID) %>% 
  mutate(fit = nls(Value ~ A * exp(-k * x), start = c(A =2, k = 0.01)) %>% list())

方法二：列表循环拟合 #

# Way2 - 构造x列表
L1 = c()
for (i in unique(Mfrq.df.2$User.ID)) {L1[[as.character(i)]]=seq(1,10)}
length(L1) # 520 users
dput(head(L1))

list(`37593` = 1:10, `38643` = 1:10, `49433` = 1:10, `60403` = 1:10, 
    `70923` = 1:10, `85363` = 1:10)

# Way2 - 构造y列表
L2=list.ids.RecSOC.2
length(L2) # 520 users
dput(head(L2))

list(`37593` = c(9L, 10L, 5L, 2L, 2L, 4L, 6L, 2L, 3L, 2L), `38643` = c(3L, 
0L, 2L, 2L, 10L, 0L, 8L, 17L, 20L, 11L), `49433` = c(4L, 29L, 
17L, 16L, 16L, 9L, 9L, 16L, 7L, 5L), `60403` = c(80L, 113L, 77L, 
47L, 40L, 22L, 8L, 24L, 30L, 11L), `70923` = c(19L, 21L, 7L, 
4L, 1L, 1L, 0L, 2L, 0L, 2L), `85363` = c(0L, 1L, 2L, 3L, 8L, 
7L, 6L, 1L, 4L, 3L))

# Way2 - 执行拟合
control=nls.control(maxiter=1000)
res <- mapply(function(x,y){
  nls(y~A*(exp(-k*x)),
      start=list(A=100, k=0.01), control=control,
      trace= TRUE, data=data.frame(x, y))},L1,L2, SIMPLIFY=FALSE)

解决方案 #

1. 为每个用户生成个性化初始值 #

奇异梯度的核心原因是初始值与最优值偏差过大，可通过线性化转换先估算初始值：

• 对指数模型 Value = A * exp(-k*x) 取对数（需保证Value>0，若有0可加极小值如1e-6）：log(Value) = log(A) -k*x
• 用lm()拟合线性模型，得到log(A)和k的初始估计，再转换为A=exp(log_A)

修改拟合函数，加入初始值自动计算，同时引入鲁棒性更强的nlsLM作为备选：

library(minpack.lm)

fit_exp_model <- function(x, y) {
  # 处理0值避免log报错
  y_safe <- ifelse(y == 0, 1e-6, y)
  # 线性化拟合求初始值
  lm_fit <- lm(log(y_safe) ~ x)
  init_A <- exp(coef(lm_fit)[[1]])
  init_k <- -coef(lm_fit)[[2]]
  # 确保k初始值为正（符合指数衰减逻辑）
  init_k <- max(init_k, 1e-4)
  
  # 优先尝试nls，失败则用nlsLM
  tryCatch({
    nls(y ~ A * exp(-k*x), start = list(A=init_A, k=init_k), 
        control = nls.control(maxiter=1000))
  }, error = function(e) {
    nlsLM(y ~ A * exp(-k*x), start = list(A=init_A, k=init_k), 
          control = nls.lm.control(maxiter=1000))
  })
}

# 批量拟合
res <- mapply(fit_exp_model, L1, L2, SIMPLIFY=FALSE)

2. 过滤不符合趋势的用户数据 #

部分用户数据可能不满足指数衰减特征（如单调递增、波动极端），可提前筛选：

# 筛选与x呈负相关的用户（符合衰减趋势）
valid_users <- names(L2)[sapply(L2, function(y) cor(1:10, y, method="spearman") < 0)]
L1_filtered <- L1[valid_users]
L2_filtered <- L2[valid_users]

# 对筛选后的数据拟合
res_filtered <- mapply(fit_exp_model, L1_filtered, L2_filtered, SIMPLIFY=FALSE)

3. 提取拟合系数 #

拟合完成后，批量提取每个用户的A和k：

coef_df <- do.call(rbind, lapply(res, function(fit) {
  coefs <- coef(fit)
  data.frame(A = coefs["A"], k = coefs["k"])
}))
coef_df$User.ID <- rownames(coef_df)

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