在R中用nls拟合误差函数至少量实验数据时遇奇异梯度矩阵问题

需要将误差函数拟合到6个实验数据点，已知实验时间t=6.504601e-05秒，仅需拟合参数D和s。初始猜测的参数绘图显示与数据吻合尚可，但使用nls拟合时始终出现奇异梯度矩阵错误。

误差函数代码 #

func_erf <- function(x, #m
                 D,     #m2/s
                 t,     #s
                 s      #m
){ 
  result = (erf((x+s)/sqrt(8*D*t)) - erf((x-s)/sqrt(8*D*t)))/(2*erf(s/sqrt(8*D*t)))
  return(result)
}

实验数据 #

> data_exp
1 1.000000000 0.000
2 0.766766619 0.001
3 0.252337795 0.002
4 0.098405369 0.003
5 0.046523446 0.004
6 0.004363998 0.005

数据的dput结构：

structure(list(y = c(1.00000000046026, 0.766766619156469, 0.252337794969704, 
0.0984053685324868, 0.0465234458835242, 0.00436399807604814), 
    x = c(0, 0.001, 0.002, 0.003, 0.004, 0.005)), class = "data.frame", row.names = c(NA, 
-6L))

拟合代码 #

coef_fit_erf_guess = c(1e-8,       #D, #m2/s
                       0.75*1e-3   #s, #m
)
t_exp = 6.504601e-05 #seconds

fit_nls_erf<- nls(y~func_erf(x,D,t= t_exp, s),data=data_exp,
                             start=list(D = coef_fit_erf_guess [1],
                                        s = coef_fit_erf_guess [2]))

疑问 #

1. 拟合流程为何失效？
2. 可以做出哪些改进？
3. 是否有方法找到更优的拟合初始值或采用手动迭代方式进行拟合？

拟合失效的原因 #

• 参数尺度差异过大：D的量级是1e-8，s是1e-3，两者尺度差达5个数量级，导致nls的梯度矩阵计算出现数值奇异，无法正常更新参数。
• 函数数值稳定性差：当D或s取某些值时，sqrt(8*D*t)可能过小或过大，使erf的自变量超出有效计算范围，梯度计算出现NaN或无穷大，进而引发奇异矩阵。
• 数据点数量不足：仅6个数据点拟合2个参数，参数自由度相对较高，算法容易陷入局部极值或数值不稳定区域。

改进措施 #

1. 参数尺度归一化 #

对参数做尺度变换，消除量级差异，将参数映射到接近1的区间后再拟合，最后反变换得到原始参数：

# 定义尺度变换后的函数
func_erf_scaled <- function(x, D_scaled, t, s_scaled) {
  D <- D_scaled * 1e-8
  s <- s_scaled * 1e-3
  result <- (erf((x+s)/sqrt(8*D*t)) - erf((x-s)/sqrt(8*D*t)))/(2*erf(s/sqrt(8*D*t)))
  return(result)
}

# 初始猜测同步变换
coef_guess_scaled <- c(1, 0.75)
fit_nls_scaled <- nls(y~func_erf_scaled(x, D_scaled, t=t_exp, s_scaled),
                      data=data_exp,
                      start=list(D_scaled=coef_guess_scaled[1], s_scaled=coef_guess_scaled[2]))

# 反变换得到原始参数
coef_fit <- coef(fit_nls_scaled)
D_fit <- coef_fit[1] * 1e-8
s_fit <- coef_fit[2] * 1e-3

2. 改用稳健拟合算法 #

放弃nls，使用minpack.lm包的nlsLM函数，它基于Levenberg-Marquardt算法，对初始值要求更低，数值稳定性更强：

library(minpack.lm)

fit_lm <- nlsLM(y~func_erf(x, D, t=t_exp, s),
                data=data_exp,
                start=list(D=1e-8, s=0.75e-3))
summary(fit_lm)

3. 优化初始值 #

• 网格搜索法：在合理范围内对D和s做网格采样，计算每个组合的残差平方和，选择残差最小的组合作为初始值：

# 定义参数搜索范围
D_grid <- seq(5e-9, 1.5e-8, by=1e-9)
s_grid <- seq(0.5e-3, 1e-3, by=0.1e-3)

# 计算所有参数组合的残差平方和
residuals_df <- expand.grid(D=D_grid, s=s_grid)
residuals_df$rss <- apply(residuals_df, 1, function(row) {
  y_pred <- func_erf(data_exp$x, row$D, t_exp, row$s)
  sum((data_exp$y - y_pred)^2)
})

# 选取最优初始值
best_guess <- residuals_df[which.min(residuals_df$rss), ]

• 物理约束缩小范围：结合实验背景，明确D和s的物理合理取值区间，避免算法进入无意义的参数区域。

4. 手动迭代拟合 #

如果自动拟合始终失败，可采用交替固定参数的手动迭代方式：

# 初始参数
current_D <- 1e-8
current_s <- 0.75e-3
rss_prev <- Inf

# 迭代次数设置
max_iter <- 10
for (i in 1:max_iter) {
  # 固定D，优化s
  fit_s <- nls(y~func_erf(x, current_D, t=t_exp, s),
               data=data_exp,
               start=list(s=current_s))
  current_s <- coef(fit_s)[["s"]]
  
  # 固定s，优化D
  fit_D <- nls(y~func_erf(x, D, t=t_exp, current_s),
               data=data_exp,
               start=list(D=current_D))
  current_D <- coef(fit_D)[["D"]]
  
  # 计算当前残差，检查收敛
  y_pred <- func_erf(data_exp$x, current_D, t_exp, current_s)
  rss_current <- sum((data_exp$y - y_pred)^2)
  
  if (abs(rss_current - rss_prev) < 1e-10) break
  rss_prev <- rss_current
  cat(sprintf("Iteration %d: D=%.2e, s=%.2e, RSS=%.4e\n", i, current_D, current_s, rss_current))
}

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