为什么来自 mgcv 的 bam 对于某些数据很慢？
Why is bam from mgcv slow for some data?
我正在使用 mgcv 中的 bam 函数在多个数据集上拟合相同的广义加性模型。而对于我的大部分数据集，拟合在 10 到 20 分钟之间的合理时间内完成。对于少数数据集，运行 需要 10 多个小时才能完成。我找不到慢案例之间的任何相似之处，最终拟合既不是特别好也不是特别差，也不包含任何明显的异常值。
我怎样才能弄清楚为什么这些实例的拟合速度如此之慢？我怎样才能加快这些速度？
我的模型包含两个平滑项（使用循环三次样条基）和一些额外的数值和因子变量。总共估计了 300 个系数（包括平滑项的系数）。我故意将结数保持在信息理论上的最佳数字以下，以加快拟合过程。我的数据集包含大约 85 万行。
这是函数调用：
bam(
value
~ 0
+ weekday_x
+ weekday
+ time
+ "a couple of factor variables encoding special events"
+ delta:weekday
+ s(share_of_year, k=length(knotsYear), bs="cc")
+ s(share_of_year_x, k=length(knotsYear), bs="cc")
, knots=list(
share_of_year=knotsYear
, share_of_year_x=knotsYear
)
, family=quasipoisson()
, data=data
)
knotsYears 包含 26 节。
这个模型在大多数情况下收敛得相当快，但在少数情况下收敛得非常慢。
最可能的原因："fREML"失败
在像上面这样的典型模型中，没有张量平滑 te 或 ti，我的经验是 REML 迭代在某些情况下会失败。
规范的 bam 实现使用 fast.REML.fit。该例程的收敛测试需要修复，但随着 Simon 继续前进并更加关注 discrete 方法，他并不热衷于修复它。固定版本（目前）仅在用于测试的扩展包中可用，"Zheyuan add-on"，作为我博士论文的补充。但是 fast.REML.fit 也不是那么脆弱，这种收敛失败的情况并不常见，否则在 2012 年一堆大报告就会解决这个问题。
以下内容只是帮助您进行检查，而不是修复。
让 fit 成为您需要 10 小时的拟合模型，检查 fit$outer.info。这给出了它需要的 REML 迭代次数，以及梯度和 Hessian 等收敛信息。如果您看到 iter = 200，或任何说明某些失败的信息，如 "step failed"，那么您就知道为什么要花这么长时间了。但是如果你看一下梯度，你很可能会发现它几乎为零。换句话说，REML 迭代实际上已经收敛，但 fast.REML.fit 未能检测到它。
另一个检查："performance iteration"
由于您拟合的是 GAM 而不是 AM（具有高斯响应的加性模型），因此在 REML 迭代之外还有另一个 P-IRLS（惩罚迭代重新加权最小二乘法）。是的，整个（规范的）bam 估计是一个双循环嵌套，称为 "performance iteration"。这也可能会失败，但这种失败更为内在，可能无法克服，因为 "performance iteration" 不能保证收敛。因此，检查 fit$iter 看它是否也很大（在最坏的情况下可以是 200）。 mgcv 手册有一个专门的部分 ?gam.convergence 讨论这种类型的收敛失败，这是 "outer iteration" 可取的驱动原因。但是，对于大型数据集，"outer iteration" 的实现成本太高。所以，忍受 "performance iteration".
mgcv 有一个 "tracing" 选项。如果在调用 bam 时设置 control = gam.control(trace = TRUE)，偏差信息和迭代计数器将在 "performance iteration" 时打印到屏幕上。这为您提供了惩罚偏差的清晰路径，因此您可以检查它是在循环还是在某个点被困。这比存储在 fit$iter.
中的单个迭代编号提供更多信息
也许试试新方法？
也许您想尝试新的 discrete = TRUE（2015 年添加；2017 年正式发表论文）。它使用新的拟合迭代。与旧方法相比，我们（更）乐于测试其实际收敛能力。使用时，也要打开"trace"。如果无法收敛，可以考虑举报，但我们需要一个可复现的案例。
我正在使用 mgcv 中的 bam 函数在多个数据集上拟合相同的广义加性模型。而对于我的大部分数据集，拟合在 10 到 20 分钟之间的合理时间内完成。对于少数数据集，运行 需要 10 多个小时才能完成。我找不到慢案例之间的任何相似之处，最终拟合既不是特别好也不是特别差，也不包含任何明显的异常值。
我怎样才能弄清楚为什么这些实例的拟合速度如此之慢？我怎样才能加快这些速度？
我的模型包含两个平滑项（使用循环三次样条基）和一些额外的数值和因子变量。总共估计了 300 个系数（包括平滑项的系数）。我故意将结数保持在信息理论上的最佳数字以下，以加快拟合过程。我的数据集包含大约 85 万行。
这是函数调用：
bam(
value
~ 0
+ weekday_x
+ weekday
+ time
+ "a couple of factor variables encoding special events"
+ delta:weekday
+ s(share_of_year, k=length(knotsYear), bs="cc")
+ s(share_of_year_x, k=length(knotsYear), bs="cc")
, knots=list(
share_of_year=knotsYear
, share_of_year_x=knotsYear
)
, family=quasipoisson()
, data=data
)
knotsYears 包含 26 节。
这个模型在大多数情况下收敛得相当快，但在少数情况下收敛得非常慢。
最可能的原因："fREML"失败
在像上面这样的典型模型中，没有张量平滑 te 或 ti，我的经验是 REML 迭代在某些情况下会失败。
规范的 bam 实现使用 fast.REML.fit。该例程的收敛测试需要修复，但随着 Simon 继续前进并更加关注 discrete 方法，他并不热衷于修复它。固定版本（目前）仅在用于测试的扩展包中可用，"Zheyuan add-on"，作为我博士论文的补充。但是 fast.REML.fit 也不是那么脆弱，这种收敛失败的情况并不常见，否则在 2012 年一堆大报告就会解决这个问题。
以下内容只是帮助您进行检查，而不是修复。
让 fit 成为您需要 10 小时的拟合模型，检查 fit$outer.info。这给出了它需要的 REML 迭代次数，以及梯度和 Hessian 等收敛信息。如果您看到 iter = 200，或任何说明某些失败的信息，如 "step failed"，那么您就知道为什么要花这么长时间了。但是如果你看一下梯度，你很可能会发现它几乎为零。换句话说，REML 迭代实际上已经收敛，但 fast.REML.fit 未能检测到它。
另一个检查："performance iteration"
由于您拟合的是 GAM 而不是 AM（具有高斯响应的加性模型），因此在 REML 迭代之外还有另一个 P-IRLS（惩罚迭代重新加权最小二乘法）。是的，整个（规范的）bam 估计是一个双循环嵌套，称为 "performance iteration"。这也可能会失败，但这种失败更为内在，可能无法克服，因为 "performance iteration" 不能保证收敛。因此，检查 fit$iter 看它是否也很大（在最坏的情况下可以是 200）。 mgcv 手册有一个专门的部分 ?gam.convergence 讨论这种类型的收敛失败，这是 "outer iteration" 可取的驱动原因。但是，对于大型数据集，"outer iteration" 的实现成本太高。所以，忍受 "performance iteration".
mgcv 有一个 "tracing" 选项。如果在调用 bam 时设置 control = gam.control(trace = TRUE)，偏差信息和迭代计数器将在 "performance iteration" 时打印到屏幕上。这为您提供了惩罚偏差的清晰路径，因此您可以检查它是在循环还是在某个点被困。这比存储在 fit$iter.
也许试试新方法？
也许您想尝试新的 discrete = TRUE（2015 年添加；2017 年正式发表论文）。它使用新的拟合迭代。与旧方法相比，我们（更）乐于测试其实际收敛能力。使用时，也要打开"trace"。如果无法收敛，可以考虑举报，但我们需要一个可复现的案例。