如何优化MCMC链中预查与即时计算因子逻辑以接近纯查表性能？

在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;
}

当n=20时，循环执行该函数1亿次耗时约2.2秒（单次波动约10ms）。由于权重因子计算耗时，预计算因子数组后得到method2，性能提升显著，但无法处理n > N_PRECOMPUTED的场景：

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;
}

相同基准下method2耗时约1.1秒，但程序需支持任意n（通常为几十），因此采用混合逻辑的method3：

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=20（理论可全查表），method3仍耗时约1.5秒，性能不及method2，需优化查表与即时计算的决策逻辑，让method3性能接近method2。

优化方案 #

1. 消除循环内分支（分支预测优化） #

循环内的条件判断会破坏CPU分支预测，导致性能损耗。通过循环拆分将查表和计算部分分开执行，让每个循环内无分支，CPU可高效流水线执行：

void method4(int n, double* x, double* y) {
    double z = 0;
    int i;
    // 先处理无法查表的部分
    const int split = n - N_PRECOMPUTED;
    const int start_split = split > 0 ? split : 0;
    for (i = 0; i < start_split; ++i) {
        const int j = n - i;
        const double factor = 1 / sqrt(j * (j + 1));
        y[i] = z - factor * j * x[i];
        z += factor * x[i];
    }
    // 再处理可查表的部分
    for (; 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;
}

2. 扩展预计算范围 #

既然n通常为几十，直接将N_PRECOMPUTED设为最大预期n+1（比如100），预计算足够大的数组。这样大部分场景下完全不需要即时计算，性能和method2一致，同时支持更大的n（超出部分再即时计算）。

3. 数学简化权重因子计算 #

观察权重因子公式：factor = 1 / sqrt(j*(j+1))，可变形为：

factor = 1 / (j * sqrt(1 + 1/j))

对于较大的j（比如j>20），用泰勒展开近似sqrt(1 + 1/j) ≈ 1 + 1/(2j) - 1/(8j²)，计算速度比直接开平方更快，且精度损失在MCMC场景下可忽略。

4. 分支预测提示 #

在条件判断前添加编译器分支预测提示（GCC/clang用__builtin_expect），告诉CPU更可能走的分支，减少流水线停顿：

const double factor = __builtin_expect(j < N_PRECOMPUTED, 1) ? FACTORS[j] : 1 / sqrt(j * (j + 1));

如果大部分场景下j都小于N_PRECOMPUTED，这个提示能大幅提升分支预测准确率。

5. 动态预计算 #

程序启动时或首次遇到较大的n时，用堆内存动态扩展FACTORS数组，将用到的j值对应的因子预计算并缓存。后续再遇到相同j直接查表，避免重复计算。

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