预先计算序列的第一部分以加速计算

Question

我必须反复评估以下函数作为 MCMC 链的一部分。它评估加权累积和并消耗大约 80% 的计算时间。

void method1(int n, double* x, double* y) {
    double z = 0;  // Weighted running total.
    for (int i = 0; i < n; ++i) {
        // Evaluate weighting factor.
        const int j = n - i;
        const double factor = 1 / sqrt(j * (j + 1));
        // Set next result.
        y[i] = z - factor * j * x[i];
        // Update running total.
        z += factor * x[i];
    }
    y[n] = z;
}

在我的机器上循环运行 100m 次

n = 20

大约需要 2.2 秒（运行之间大约有 10 毫秒的变化）。

这里“昂贵”的操作是计算权重因子，我可以在启动 MCMC 链之前预先计算因子。

double FACTORS[N_PRECOMPUTED] = ...;

void method2(int n, double* x, double* y) {
    double z = 0;
    for (int i = 0; i < n; ++i) {
        const int j = n - i;
        const double factor = FACTORS[j];
        y[i] = z - factor * j * x[i];
        z += factor * x[i];
    }
    y[n] = z;
}

运行相同的基准测试，大约需要 1.1 秒，这似乎是一个值得的改进。然而，当

n > N_PRECOMPUTED

时，这当然会爆炸。虽然

n

通常很小（十的数量级），但程序应该针对任何

n

运行。

下一个逻辑步骤是预先计算序列的第一部分，并根据需要动态评估其余部分。

void method3(int n, double* x, double* y) {
    double z = 0;
    for (int i = 0; i < n; ++i) {
        const int j = n - i;
        const double factor = j < N_PRECOMPUTED ? FACTORS[j] : evaluate_factor(j);
        y[i] = z - factor * j * x[i];
        z += factor * x[i];
    }
    y[n] = z;
}

这对

N_PRECOMPUTED = 10

产生了一些改进：1.8 秒，但不如纯查找减少 50% 那样令人印象深刻。如果我设置

N_PRECOMPUTED = 20

，原则上它可以查找所有值，它会在 1.5 秒内运行。

是否有更好的方法来处理在

method3

中动态查找或计算因子的决策，以更接近

method2

的性能？

Answer 1

也许可以使用混合策略，在其中动态计算一些因素，同时尽可能地利用预先计算的因素。这是方法3的修改版本：

void method4(int n, double* x, double* y) {
    double z = 0;
    int i = 0;
    
    // Handle the initial precomputed factors if available
    for (; i < N_PRECOMPUTED && i < n; ++i) {
        const int j = n - i;
        const double factor = FACTORS[j];
        y[i] = z - factor * j * x[i];
        z += factor * x[i];
    }

    // Calculate and store remaining factors on-the-fly
    for (; i < n; ++i) {
        const int j = n - i;
        const double factor = evaluate_factor(j);
        y[i] = z - factor * j * x[i];
        z += factor * x[i];
    }

    y[n] = z;
}

此方法首先迭代预先计算的因子（如果可用），然后即时计算剩余因子。这样，您就可以在可能的情况下利用预先计算的因素，并在必要时回退到即时计算。确切的性能改进将取决于n值的分布以及可以预先计算多少个因素，但它应该比方法3更有效。

您可以根据典型用例和可用内存调整 N_PRECOMPUTED，以找到预计算和即时计算之间的最佳权衡。