CRAN 包提交：“错误：C 堆栈使用量太接近限制”

Question

正确的：这是我在向 CRAN 提交 R 包时遇到的问题。所以我

• 无法控制堆栈大小（因为问题发生在某个 CRAN 平台上）
• 我无法提供可重现的示例（因为我不知道 CRAN 上的确切配置）

问题

当尝试将cSEM.DGP package 提交给 CRAN 时，自动预测试（适用于 Debian x86_64-pc-linux-gnu；不适用于 Windows！）失败并出现注释：C stack usage 7975520 is too close to the limit。

我知道这是由一个具有三个参数的函数引起的，该函数的主体大约有 800 行长。函数体由这些参数的加法和乘法组成。它是函数 varzeta6()，您可以找到 here（从第 647 行开始）。

我该如何解决这个问题？

我不能做的事情：

• 提供一个可重现的例子（至少我不知道怎么做）
• 更改堆栈大小

我在想的事情：

• 尝试将函数分解为更小的部分。但我不知道如何最好地做到这一点。
• 以某种方式预编译？功能（老实说，我只是猜测）所以 CRAN 不会抱怨？

让我知道你的想法！

细节/背景

varzeta6()（以及varzeta4()/varzeta5() 甚至更多varzeta7()）如此长且 R 效率低的原因是它们本质上是从mathematica 复制粘贴的（在将mathematica 代码简化为尽可能好并将其调整为有效的 R 代码）。因此，代码绝不是 R 优化的（@MauritsEvers 正确地指出了这一点）。

为什么我们需要mathematica？因为我们需要的是递归结构方程模型的模型隐含构造相关矩阵的一般形式，该模型具有多达 8 个构造作为模型方程参数的函数。此外还有约束。 为了了解这个问题，让我们采用一个可以递归求解的两个方程组：

• Y2 = beta1*Y1 + zeta1
• Y3 = beta2*Y1 + beta3*Y2 + zeta2

我们感兴趣的是协方差：E(Y1*Y2)、E(Y1*Y3) 和 E(Y2*Y3) 作为 beta1、beta2、beta3 的函数，在

• E(Y1) = E(Y2) = E(Y3) = 0,
• E(Y1^2) = E(Y2^2) = E(Y3^3) = 1
• E(Yi*zeta_j) = 0（其中 i = 1、2、3 和 j = 1、2）

对于这样一个简单的模型，这是相当微不足道的：

• E(Y1*Y2) = E(Y1*(beta1*Y1 + zeta1) = beta1*E(Y1^2) + E(Y1*zeta1) = beta1
• E(Y1*Y3) = E(Y1*(beta2*Y1 + beta3*(beta1*Y1 + zeta1) + zeta2) = beta2 + beta3*beta1
• E(Y2*Y3) = ...

但是，当您添加 Y4、Y5 直到 Y8 时，您会发现这很快就会变得混乱。 一般来说，模型隐含的构造相关矩阵可以写成（表达式实际上看起来更复杂，因为我们还允许多达 5 个外生构造。这就是为什么varzeta1() 已经看起来很复杂。但现在忽略这个。） ：

• V(Y) = (I - B)^-1 V(zeta)(I - B)'^-1

其中 I 是单位矩阵，B 是模型参数（β）的下三角矩阵。 V(zeta) 是对角矩阵。函数varzeta1()、varzeta2()、...、varzeta7() 计算主对角线元素。由于我们将 Var(Yi) 限制为始终为 1，因此 zeta 的方差随之而来。以方程 Var(Y2) = beta1^2*Var(Y1) + Var(zeta1) --> Var(zeta1) = 1 - beta1^2 为例。这在这里看起来很简单，但是当我们取这样一个递归方程链中的第 6 个方程的方差时，它变得极其复杂，因为 Var(zeta6) 取决于所有先前的 Y1、...、Y5 之间的协方差，它们本身就是取决于它们各自先前的协方差。

好吧，我不知道这是否能让事情变得更清楚。以下是重点：

1. varzeta1(), ..., varzeta7() 的代码是从 Mathematica 复制粘贴的，因此未经过 R 优化。
2. Mathematica 是必需的，因为据我所知，R 无法处理符号计算。
3. 我可以“手动”进行 R 优化（这非常乏味）
4. 我认为varzetaX() 的结构必须是给定的。因此，问题是：我可以以某种方式使用这个功能吗？

Answer 1

这是一个有点太长的评论，但希望这会给你一些优化varzeta*函数代码的想法；或者至少，它可能会给你一些思考。

有几件事让我感到困惑：

1. 所有varzeta* 函数都有参数beta、gamma 和phi，它们似乎是矩阵。但是，在varzeta1 中，您不使用beta，但beta 是第一个函数参数。
2. 我很难将您在帖子底部提供的详细信息与varzeta* 函数的代码联系起来。您没有解释 gamma 和 phi 矩阵的来源，也没有解释它们的含义。此外，看到beta 是模型的参数估计，我不明白为什么beta 应该是一个矩阵。

正如我在之前的评论中提到的，如果这些表达式不能简化，我会感到非常惊讶。 R 可以很轻松地进行很多矩阵运算，实际上不需要预先计算各个项。

例如可以使用crossprod和tcrossprod计算叉积，%*%实现矩阵乘法。

其次，R 中的许多数学运算都是矢量化的。我已经提到你可以简化

1 - gamma[1,1]^2 - gamma[1,2]^2 - gamma[1,3]^2 - gamma[1,4]^2 - gamma[1,5]^2

作为

1 - sum(gamma[1, ]^2)

因为^ 运算符是矢量化的。

---

也许从根本上说，这对我来说似乎有点 XY 问题，这可能有助于退后一步。不知道您要建模的全部细节（正如我所说，我无法将您提供的细节链接到 cSEM.DGP 代码），我将首先探索如何解决 R 中的递归 SEM。在这里并没有真正看到对 Mathematica 的需求。正如我之前所说，矩阵运算在 R 中是非常标准的；在 R 中也可以解析求解一组递归方程。由于您似乎来自 Mathematica 领域，因此最好与当地的 R 编码专家讨论这个问题。

如果你必须使用那些可怕的 varzeta* 函数（我真的怀疑），一个选项可能是用 C++ 重写它们，然后用 Rcpp 编译它们以将它们变成 R职能。也许这样可以避免 C 堆栈使用限制？

Answer 2

一种可行的方法是尝试说服 CRAN 维护人员，让您相信没有简单的方法可以解决问题。这是NOTE，而不是WARNING； CRAN repository policy 说

原则上，包裹必须在没有警告或重要说明的情况下通过 R CMD 检查才能进入主 CRAN 包裹区域。如果存在您无法消除的警告或注释（例如，因为您认为它们是虚假的），请发送解释性说明作为您的覆盖电子邮件的一部分，或作为对提交表单的评论

因此，您可以把握机会，您的合理解释（在提交表单上的 cmets 字段中）将说服 CRAN 维护人员。从长远来看，最好找到一种方法来简化计算，但在提交给 CRAN 之前可能没有必要这样做。