Implementing the Shekel function in Matlab

Question

I am trying to write a faster way to evaluate the Shekel function found here : https://www.sfu.ca/~ssurjano/shekel.html

Their code is pretty similar to the one I have been using, except I pass a matrix of x values to the function and end up with three loops in a element by element calculation. Matlab should be able to do better.

This is the code I had been using previously :

function S = Shekel(X,m)
[R,d] = size(X);
% R is the population size, m the number of minima, and d the dimensions

% input control %
if d > 4
    error('More than 4 dimensions !!')
end
if nargin==1
    m=10;
elseif (m > 10) || (m < 2)
    error('Wrong m')
end

% coefficients %
A = [4 4 4 4
     1 1 1 1
     8 8 8 8
     6 6 6 6
     3 7 3 7
     2 9 2 9
     5 5 3 3
     8 1 8 1
     6 2 6 2
     7 3.6 7 3.6];

c = [.1 .2 .2 .4 .4 .6 .3 .7 .5 .5];

S = zeros(R,1);
for r = 1:R
    s = 0;
    for i = 1:m
        denom = c(i);
        for j = 1:d
            denom = denom + (X(r,j) - A(i,j))^2;
        end
        s = s - 1/denom;
    end
    S(r) = s;
end

There is also this R implementation which is similar to what I would like to do : https://www.sfu.ca/~ssurjano/Code/shekelr.html

So far I have gotten to this, but it is as slow as the previous one because I'm doing exactly the same thing :

S = zeros(R,1);
for r = 1:R
    %S(r) = -sum(1./(sum((repmat(X(r,:),m,1)-A).^2,2)'+c),2);
    S(r) = -sum(1./(sum((repmat(X(r,:),m,1)-A).^2,2).'+c),2);
end

The function can be called with this and you should get -0.3007 :

clear

p = 40; m = 10; d = 4;
X = repmat([1 2 3 4], p, 1);
OF = @Shekel;
OF(X, m)

I call Shekel a lot and it is responsible for most of my execution time. Any suggestions ?

_{EDIT: Actual operations (population size = 2, dimensions = 4, #minimums = 10)}

_{X =}

1.5381    0.7603    3.2619    7.7624
8.1874    1.9172    0.4234    5.0153

_{A =}

4.0    4.0    4.0    4.0
1.0    1.0    1.0    1.0
8.0    8.0    8.0    8.0
6.0    6.0    6.0    6.0
3.0    7.0    3.0    7.0
2.0    9.0    2.0    9.0
5.0    5.0    3.0    3.0
8.0    1.0    8.0    1.0
6.0    2.0    6.0    2.0
7.0    3.6    7.0    3.6

_{repmat(X(1,:),m,1)-A}

1.5381    0.7603    3.2619    7.7624        4.0    4.0    4.0    4.0
1.5381    0.7603    3.2619    7.7624        1.0    1.0    1.0    1.0
1.5381    0.7603    3.2619    7.7624        8.0    8.0    8.0    8.0
1.5381    0.7603    3.2619    7.7624        6.0    6.0    6.0    6.0
1.5381    0.7603    3.2619    7.7624    -   3.0    7.0    3.0    7.0
1.5381    0.7603    3.2619    7.7624        2.0    9.0    2.0    9.0
1.5381    0.7603    3.2619    7.7624        5.0    5.0    3.0    3.0
1.5381    0.7603    3.2619    7.7624        8.0    1.0    8.0    1.0
1.5381    0.7603    3.2619    7.7624        6.0    2.0    6.0    2.0
1.5381    0.7603    3.2619    7.7624        7.0    3.6    7.0    3.6

_{X-A =}

-2.4619   -3.2397   -0.7381    3.7624
 0.5381   -0.2397    2.2619    6.7624
-6.4619   -7.2397   -4.7381   -0.2376
-4.4619   -5.2397   -2.7381    1.7624
-1.4619   -6.2397    0.2619    0.7624
-0.4619   -8.2397    1.2619   -1.2376
-3.4619   -4.2397    0.2619    4.7624
-6.4619   -0.2397   -4.7381    6.7624
-4.4619   -1.2397   -2.7381    5.7624
-5.4619   -2.8397   -3.7381    4.1624

_{(repmat(X(1,:),m,1)-A).^2

(X-A)^2}

 6.0612   10.4955     0.5448    14.1560
 0.2895    0.0574     5.1162    45.7307
41.7568   52.4129    22.4495     0.0564
19.9090   27.4542     7.4971     3.1062
 2.1373   38.9335     0.0686     0.5813
 0.2134   67.8922     1.5924     1.5315
11.9851   17.9748     0.0686    22.6809
41.7568    0.0574    22.4495    45.7307
19.9090    1.5368     7.4971    33.2058
29.8329    8.0638    13.9733    17.3260

_{sum((repmat(X(1,:),m,1)-A).^2,2).'

sum of X-A rows transposed}

31.2575  51.1939  116.6756  57.9665  41.7208  71.2296  52.7094  109.9944  62.1487  69.1959

_{c =}

0.1  0.2  0.2  0.4  0.4  0.6  0.3  0.7  0.5  0.5

_{sum((repmat(X(1,:),m,1)-A).^2,2).'+c

sum(X-A)+c}

31.3575  51.3939  116.8756  58.3665  42.1208  71.8296  53.0094  110.6944  62.6487  69.6959

_{1./(sum((repmat(X(r,:),m,1)-A).^2,2).'+c)

1/sum(X-A)+c}

0.0319  0.0195  0.0086  0.0171  0.0237  0.0139  0.0189  0.0090  0.0160  0.0143

_{-sum(1./(sum((repmat(X(1,:),m,1)-A).^2,2).'+c),2)

Shekel(X(1)) = -sum(1/sum(X-A)+c)}

-0.1729

_{.. and so on for X(2)}

You can probably remove the `repmat` for a small speed gain. `-` will do *singleton expansion`, meaning that `repmat` will be implicitly applied to match up the sizes of the two operands. Because it's implicit it uses less memory, and hence the computations are quicker. — Cris Luengo, Feb 25 '19 at 16:02
I changed to `.'` (transpose), but `X(r,:)-A` without `repmat()` isn't working `(Matrix dimensions must agree.)` and `S(r) = -sum(1./(sum((bsxfun(@minus,X(r,:),A)).^2,2).'+c),2);` for some reasons is 2x slower. Maybe I should add an exploded view on what operations are being performed to see if there's a more 'elegant' way or if elementwise is forced. — alfred, Feb 26 '19 at 13:54
If removing `repmat` doesn’t work, you’re using an older version of MATLAB. The implicit singleton expansion works in R2016b and newer I think. — Cris Luengo, Feb 26 '19 at 14:13

Implementing the Shekel function in Matlab

0 Answers0