Optimizing MATLAB code

Question

This code takes an extremely long time to run (more than 10 minutes). Is there any way in which I can optimize it so that it finishes in less than one minute?

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

The inner loop is executing around 1000000*(1000000+1)/2 = 500000500000 times! No wonder it is slow. Maybe you should try a different approximation approach.

Answer 2

@b3 has a great vectorization of rhs.

One typo though, needs to use times and not mtimes:

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

For lhs, I propose the following, loosely based on Eratosthenes' Sieve:

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;

Execution time is just 2.45 seconds (for this half of the problem). Total including calculation of rhs and find is under 3 seconds.

I'm currently running the other version to make sure that the results are equal.

---

EDIT: found a bug with lhs(1) and special-cased it (it is a special case, the only natural number where 1 and N aren't distinct factors)

Answer 3

Vectorizing your algorithm where I could reduced the execution time slightly to ~8.5 minutes. Calculate all of the harmonic sums in one statement:

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

You can now calculate the right-hand side in one statement:

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

I couldn't vectorize the determination of the factors so this is the fastest way I could come up with to sum them. MATLAB's FACTOR command allows you to generate all the prime factors for each iteration. We then compute the unique set of products of all possible combinations using UNIQUE and NCHOOSEK. This avoids testing each integer as a factor.

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;

Find the indices at which the Riemann hypothesis is violated:

isViolated = find(lhs > rhs);