Оптимизация кода 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

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 не являются различными факторами)

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

Answer 3

Внутренний цикл выполняется около 1000000 * (1000000 + 1) / 2 = 500000500000 раз! Не удивительно, что это медленно. Может быть, вам следует попробовать другой подход приближения.