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, Oct 03 '11 at 22:06

score 6 · Accepted Answer · edited May 23 '17 at 12:19

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 '17 at 12:19

Community

1
1

answered Oct 03 '11 at 21:29

Ben Voigt

277,958
43
419
720

Thanks for catching the typo. Nice job on the `lhs` calculation. I think all you need is to add 1 to the entire vector. – b3. Oct 03 '11 at 21:54
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? – icobes Oct 03 '11 at 22:34
@icobes: It's the magic of `cumsum`. Note that `harmonicsum(i) = harmonicsum(i-1) + 1/i`. – Ben Voigt Oct 03 '11 at 22:36
+1 this is seriously fast.. it computes the [sum of divisors](http://oeis.org/A000203) 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... – Amro Oct 04 '11 at 02:59

b3. · Answer 2 · 2011-10-03T22:24:09.610

4

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

edited Oct 03 '11 at 22:24

answered Oct 03 '11 at 21:05

b3.

7,094
2
33
48

As far as I can tell, you get wrong results for `lhs(1)`, `lhs(2)`, `lhs(3)`, and probably all the rest also. – Ben Voigt Oct 03 '11 at 21:58
Ran the original code from the question to find `lhs(1:10)`, mine now matches, yours doesn't, not for a single element :( – Ben Voigt Oct 03 '11 at 22:03
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. – Ben Voigt Oct 03 '11 at 22:08
That should have said "change the loop to `for jj = 1:numFactor-1`". – Ben Voigt Oct 03 '11 at 22:16
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? – icobes Oct 03 '11 at 22:18
@icobes - Look at Ben Voigt's solution. It's much, much faster. – b3. Oct 03 '11 at 22:22
@Ben Voigt - Fixed the loop count although it is unimportant considering the faster solution. – b3. Oct 03 '11 at 22:25

score 0 · Answer 3 · edited Dec 27 '11 at 21:21

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 '11 at 21:21

Peter Mortensen

30,738
21
105
131

answered Oct 03 '11 at 20:08

Óscar López

232,561
37
312
386

I don't know how else to test for factors and test each of the 1,000,000 numbers without using a double loop. – icobes Oct 03 '11 at 20:21

Optimizing MATLAB code

3 Answers3

Linked