优化MATLAB代码

Question

这段代码运行时间非常长（超过10分钟），有没有什么方法可以进行优化，使得它在不到一分钟内完成？

clear all;
for i = 1:1000000
    harmonicsum = 0;
    lhs = 0;
    for j = 1:i
        % compute harmonic sum
        harmonicsum = harmonicsum + 1/j;
        % find sum of factors
        if (mod(i,j)==0)
            lhs = lhs + j;
        end
    end
    %define right hand side (rhs) of Riemann Hypothesis
    rhs = harmonicsum + log(harmonicsum) * exp(harmonicsum);

    if lhs > rhs
        disp('Hypothesis violated')
    end
end

- icobes

---

关键不是问“N的因数是什么”，而是“哪些数字有j作为因数”，这要容易得多。 - Ben Voigt

Answer 1

@b3有一个很好的向量化rhs的方法。

不过，有一个错别字，需要使用times而不是mtimes：

harmonicsum = cumsum(1 ./ (1:1e6));
rhs = harmonicsum + log(harmonicsum) .* exp(harmonicsum);

对于lhs，我建议以下内容，基于埃拉托斯特尼筛法的基本原理：

lhs = 1 + [1:1e6];
lhs(1) = 1;
for iii = 2:numel(lhs)/2
    lhs(2*iii:iii:end) = lhs(2*iii:iii:end) + iii;
end;

执行时间仅为2.45秒（针对该问题的一半）。总计包括rhs和find的计算在内不到3秒。

我目前正在运行另一个版本，以确保结果相等。

---

编辑：发现了一个与lhs（1）有关的漏洞，并对其进行了特殊处理（这是一种特殊情况，唯一的自然数，在该情况下1和N不是不同的因数）

- Ben Voigt

---

感谢您发现了这个笔误。lhs计算得很好。我认为您只需要将整个向量加1即可。 - b3.

对于谐和级数，我需要将其放在循环内吗？我不明白你的代码如何为每个i的值进行谐和求和？ - icobes

@icobes：这就是cumsum的魔力。请注意，harmonicsum(i) = harmonicsum(i-1) + 1/i。 - Ben Voigt

+1 这真的很快...它计算了n的所有因子之和（http://oeis.org/A000203），在1到1000000之间的所有n中，仅用了不到3秒钟！我尝试了该页面上列出的一些代码，Mathematica版本需要15秒，MuPad版本需要约100秒... - Amro

Answer 2

将您的算法向量化后，执行时间略微缩短至约8.5分钟。在一个语句中计算所有谐和数：

harmonicsum = cumsum(1 ./ (1:1e6));

现在，您可以在一个语句中计算右侧：

rhs = harmonicsum + log(harmonicsum) .* exp(harmonicsum);

我无法将因子确定向量化，所以这是我能想到的最快的求和方法。MATLAB的 FACTOR 命令允许您为每次迭代生成所有质因数。然后，我们使用 UNIQUE 和 NCHOOSEK 计算所有可能组合的唯一乘积集。这避免了测试每个整数作为因子。

lhs = zeros(1e6, 1);
for ii = 1:1e6
    primeFactor = factor(ii); // 获取质因数
    numFactor = length(primeFactor); // 计算质因数个数
    allFactor = [];
    for jj = 1:numFactor-1
       allFactor = [allFactor; unique(prod(nchoosek(primeFactor, jj), 2))]; // 计算所有可能的因子组合并去重
    end
    lhs(ii) = sum(allFactor) + 1 + ii; // 计算左式
end
lhs(1) = 1;

查找瑞利曼假设被违反的索引：

isViolated = find(lhs > rhs);

- b3.

---

据我所知，您对lhs（1），lhs（2），lhs（3）和可能的所有其他结果都得到了错误的结果。 - Ben Voigt

运行了问题中的原始代码以查找 lhs(1:10)，我的结果现在匹配了，而你的却没有，甚至一个元素都不匹配 :( - Ben Voigt

我认为问题（除了ii == 1的情况）在于jj == numFactor会得到ii本身，而你又单独计算了它。将循环改为jj = 1:numFactor，应该会更好。 - Ben Voigt

我的老师能在不到一分钟的时间内运行这段代码。你能想到其他表达这段代码的方式吗？ - icobes

@icobes - 看看Ben Voigt的解决方案。它快得多，快得多。 - b3.

Answer 3

内部循环执行了大约1000000 *（1000000 + 1）/ 2 = 500000500000次！难怪它很慢。也许你应该尝试不同的近似方法。

- Óscar López

---

我不知道如何在不使用双重循环的情况下测试因子并测试这100万个数字中的每一个。 - icobes