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

The key is not to ask "what are the factors of N", but "what numbers have j as a factor", which is MUCH easier. — Ben Voigt
– Ben Voigt, Commented Oct 3, 2011 at 22:06

Community · Accepted Answer · 2017-05-23 12:19:56Z

6

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

edited May 23, 2017 at 12:19

CommunityBot

11 silver badge

answered Oct 3, 2011 at 21:29

Ben Voigt

286k45 gold badges443 silver badges763 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

b3. Over a year ago

Thanks for catching the typo. Nice job on the lhs calculation. I think all you need is to add 1 to the entire vector.

icobes Over a year ago

For the harmonic sum, do I need to put it within the loop? I don't follow how your code does a harmonic sum for each value of i?

Ben Voigt Over a year ago

@icobes: It's the magic of cumsum. Note that harmonicsum(i) = harmonicsum(i-1) + 1/i.

Amro Over a year ago

+1 this is seriously fast.. it computes the sum of divisors of n, for all n between 1 and 1000000 in under 3 seconds! I tried some of the codes listed on that page, Mathematica version took 15 sec, MuPad version around 100 seconds...

b3. · Accepted Answer · 2011-10-03 22:24:09Z

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

As far as I can tell, you get wrong results for lhs(1), lhs(2), lhs(3), and probably all the rest also.
Ran the original code from the question to find lhs(1:10), mine now matches, yours doesn't, not for a single element :(
I think the problem (apart from the ii == 1 case) is that jj == numFactor yields ii itself, which you're counting separately. Change the loop to jj = 1:numFactor and it should be better.
That should have said "change the loop to for jj = 1:numFactor-1".
My teacher was able to run this code in less than one minute. Is there any other way to express this code that you can think of?

Peter Mortensen · Accepted Answer · 2011-12-27 21:21:40Z

0

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.

edited Dec 27, 2011 at 21:21

Peter Mortensen

31.4k22 gold badges110 silver badges134 bronze badges

answered Oct 3, 2011 at 20:08

Óscar López

237k38 gold badges321 silver badges391 bronze badges

1 Comment

icobes Over a year ago

I don't know how else to test for factors and test each of the 1,000,000 numbers without using a double loop.

Collectives™ on Stack Overflow

Optimizing MATLAB code

3 Answers 3

4 Comments

7 Comments

1 Comment

Linked

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

4 Comments

7 Comments

1 Comment

Linked

Related