Cartesian product of two vectors in Julia
Asked Answered
I

4

13

I have two vectors x and y, of respective lengths n and p. Is there a built-in way to create a np x 2 matrix which would be

x[1] y[1]
x[1] y[2]
...
x[1] y[p]
x[2] y[1]
...
x[n] y[p]

I can do that with a nested for loop, but I'm looking for a built-in function, if it exists.

Indeterminate answered 30/3, 2015 at 13:31 Comment(1)
Does anyone know how to modify the answer would be if I wanted by a n by p matrix of ordered pairs, where entry i j is (x[i], y[j])Fund
W
10

Julia is usually very fast in nested loops, so if the they are working correctly for you, you should proabably check the performance maybe just stick with it.

Other option would be using repmat (this one is a little faster than using repeat):

[repmat(x,1,length(y))'[:] repmat(y,length(x),1)[:]]

Did some quick testing of both methods:

x=rand(1000)
y=rand(1000)

function withrepeat(x,y)
    [repeat(x, inner=[size(y,1)]) repeat(y, outer=[size(x,1)])]
end

function withrepmat(x,y)
    [repmat(x,1,length(y))'[:] repmat(y,length(x),1)[:]]
end

withrepeat(x,y)
elapsed time: 0.21556302 seconds (95986112 bytes allocated)

with repmat(x,y)
elapsed time: 0.075604488 seconds (56000560 bytes allocated)

Not sure why so much difference and I think there is room for improvement still. Haven't tried the product function inside the Iterators.jl package.

Also a little bit more information here: https://groups.google.com/forum/#!topic/julia-users/dtl--SyFgwY

Hope this helps.

Tried a couple of nested loops and indeed is faster:

function withloops (x,y)
    leny=length(y)
    lenx=length(x)
    m=leny*lenx
    OUT = zeros(Float64, m,2)
    c=1
    for i = 1:lenx
        for j = 1:leny
            OUT[c,1] = x[i]
            OUT[c,2] = y[j]
            c+=1
        end
    end 
    return OUT
end

And, for the same rand(1000) for x and y.

withloops(x,y)
elapsed time: 0.011350679 seconds (16000128 bytes allocated)
Windle answered 30/3, 2015 at 18:38 Comment(3)
Does @inbounds help the for loop?Appliance
@Appliance It should yes, I was in a little rush when I posted it, so I didn't have time to test it properly, but I will do it tomorrow and edit the post accordingly. Thanks for the suggestion.Windle
I tried using @inbounds but didn't get any noticeable improvement, maybe this is because there are no operations with the arrays and I'm just assigning values? Also, in my computer using vectors of length 10000 as input for x and y, produced memory errors when using repeat or repmat but worked fine with the nested loops function.Windle
V
19

At least in Julia 1.3, there is a built-in function Iterators.product, so Iterators.product(a, b) |> collect should do the trick.

Vincenzovincible answered 21/11, 2019 at 9:29 Comment(1)
Note that Iterators is part of Julia's Base module (Base.Iterators), so you don't have to call import or using at all to be able to use Iterators.product.Tracheotomy
B
12

This is provided in the IterTools module (which replaces Iterators module).

Taken from https://github.com/JuliaCollections/IterTools.jl

Iterate over all combinations in the cartesian product of the inputs.

Example:

using IterTools
for p in product(1:3,1:2)
    @show p
end

p = (1,1)
p = (2,1)
p = (3,1)
p = (1,2)
p = (2,2)
p = (3,2)
Bergwall answered 9/9, 2016 at 9:19 Comment(2)
I'm not sure it replaces the Iterators module, because when using IterTools.product in julia 1.3rc4, it says it is deprecated in favour of Iterators.productVarini
This solution worked much faster for my using Julia 1.5Uro
W
10

Julia is usually very fast in nested loops, so if the they are working correctly for you, you should proabably check the performance maybe just stick with it.

Other option would be using repmat (this one is a little faster than using repeat):

[repmat(x,1,length(y))'[:] repmat(y,length(x),1)[:]]

Did some quick testing of both methods:

x=rand(1000)
y=rand(1000)

function withrepeat(x,y)
    [repeat(x, inner=[size(y,1)]) repeat(y, outer=[size(x,1)])]
end

function withrepmat(x,y)
    [repmat(x,1,length(y))'[:] repmat(y,length(x),1)[:]]
end

withrepeat(x,y)
elapsed time: 0.21556302 seconds (95986112 bytes allocated)

with repmat(x,y)
elapsed time: 0.075604488 seconds (56000560 bytes allocated)

Not sure why so much difference and I think there is room for improvement still. Haven't tried the product function inside the Iterators.jl package.

Also a little bit more information here: https://groups.google.com/forum/#!topic/julia-users/dtl--SyFgwY

Hope this helps.

Tried a couple of nested loops and indeed is faster:

function withloops (x,y)
    leny=length(y)
    lenx=length(x)
    m=leny*lenx
    OUT = zeros(Float64, m,2)
    c=1
    for i = 1:lenx
        for j = 1:leny
            OUT[c,1] = x[i]
            OUT[c,2] = y[j]
            c+=1
        end
    end 
    return OUT
end

And, for the same rand(1000) for x and y.

withloops(x,y)
elapsed time: 0.011350679 seconds (16000128 bytes allocated)
Windle answered 30/3, 2015 at 18:38 Comment(3)
Does @inbounds help the for loop?Appliance
@Appliance It should yes, I was in a little rush when I posted it, so I didn't have time to test it properly, but I will do it tomorrow and edit the post accordingly. Thanks for the suggestion.Windle
I tried using @inbounds but didn't get any noticeable improvement, maybe this is because there are no operations with the arrays and I'm just assigning values? Also, in my computer using vectors of length 10000 as input for x and y, produced memory errors when using repeat or repmat but worked fine with the nested loops function.Windle
A
6

Here's how I might do it:

julia> x = [1, 2, 3, 4];

julia> y = [9, 8, 7];

julia> [repeat(x, inner=[size(y,1)]) repeat(y, outer=[size(x,1)])]
12x2 Array{Int64,2}:
 1  9
 1  8
 1  7
 2  9
 2  8
 2  7
 3  9
 3  8
 3  7
 4  9
 4  8
 4  7

you may also want to take a look at Iterators.j -- specifically the product function.

Anglin answered 30/3, 2015 at 14:16 Comment(1)
Note that you can omit the brackets around the inner and outer arguments, so the formula reduces to [repeat(x, inner=size(y,1)) repeat(y, outer=size(x,1))]Atom

© 2022 - 2024 — McMap. All rights reserved.