Hello there,

I was wondering if you knew how to efficiently compute the following matrices in Julia. In Matlab I am getting the result 3-6x faster than in Julia, which makes me think that there is room for improvement perhaps? Since I need to compute these thousands of times in my code, I am really interested to optimize this bit…

So, here is the computation in Matlab. I take advantage of the fact that I can perform the “all-against-all” summation in matlab by summing a Nx1 with a 1xN vector, which gives the NxN matrix. In Julia, as far as I know, it is not possible.

```
function h = myfunc(v,vp)
% Both v and vp are column vectors
f = (v.^(1/3)+(vp.^(1/3))').^2.*((v.^(2/9)+(vp.^(2/9))').^(1/2));
g = exp(-(((v*vp').^(1/3))./(v.^(1/3)+(vp.^(1/3))')).^4);
h = f.*g;
end
```

I test the function with the following:

```
x=linspace(1,1000,1000)';
tic
[f, g] = myfun(x,x);
toc
```

Result:

```
Elapsed time is 0.093312 seconds.
```

Now, in Julia, since I cannot perform the NxN summation this way, I can only think of a for-loop. So I use two functions:

```
function product(v,vp)
f = zeros(Float64,length(v),length(vp));
g = zeros(Float64,length(v),length(vp));
for i in eachindex(v), j in eachindex(vp)
@inbounds f[j,i] = (v[i]^(1/3)+(vp[j]^(1/3)))^2 *((v[i]^(2/9)+(vp[j]^(2/9)))^(1/2))*2^(1/3);
@inbounds g[j,i] = exp(-(((v[i]*vp[j])^(1/3))/(v[i]^(1/3)+(vp[j]^(1/3))))^4);
end
return f, g
end
```

```
function myfun(v,vp)
f, g = product(v,vp);
h = @. f*g;
return h
end
```

I test with the same input I used in Matlab (`x=collect(1:1000)*1.0`

) and, using BenchmarkTools I get:

```
596.179 ms (8 allocations: 22.89 MiB)
```

So Matlab is 6.3 times faster.

Now, I realize that, since I am evaluating this function using the same vector for both inputs, the result is a symmetric matrix. So I can improve things using LinearAlgebra:

```
function product2(v,vp)
f = zeros(Float64,length(v),length(vp));
g = zeros(Float64,length(v),length(vp));
N = 1;
for i in eachindex(v)
for j = 1:N
@inbounds f[j,i] = (v[i]^(1/3)+(vp[j]^(1/3)))^2 *((v[i]^(2/9)+(vp[j]^(2/9)))^(1/2))*2^(1/3);
@inbounds g[j,i] = exp(-(((v[i]*vp[j])^(1/3))/(v[i]^(1/3)+(vp[j]^(1/3))))^4);
end
N = N+1;
end
f = Symmetric(f);
g = Symmetric(g);
return f, g
end
```

which gives:

```
293.099 ms (10 allocations: 22.89 MiB)
```

Better, but still 3x slower than Matlab.

I would really appreciate any thoughts on this.