# Outer product of more than two vectors

Is there a built-in function for an outer product of more than two vectors, i.e.:

``````julia> a  = [0.2, 0.8];

julia> b  = [0.1,0.9];

julia> c  = [0.2,0.2,0.6];

julia> vs = [a,b,c]; # vector of vectors

julia> function outer_product(vs) # vs is a vector of vectors
nd     = length(vs)
asize  = length.(vs)
out    = Array{Float64,nd}(undef,asize...)
for idx in CartesianIndices((asize...,))
out[idx] = *([vs[d][idx[d]] for d in 1:nd]...)
end
return out
end
outer_product (generic function with 1 method)

julia> outer_product(vs)
2×2×3 Array{Float64, 3}:
[:, :, 1] =
0.004  0.036
0.016  0.144

[:, :, 2] =
0.004  0.036
0.016  0.144

[:, :, 3] =
0.012  0.108
0.048  0.432
``````
``````julia> (a * b') .* reshape(c, 1, 1, :)
2×2×3 Array{Float64, 3}:
[:, :, 1] =
0.004  0.036
0.016  0.144

[:, :, 2] =
0.004  0.036
0.016  0.144

[:, :, 3] =
0.012  0.108
0.048  0.432
``````

Thank you, but when I don’t know the number of vectors ?

you just have to reshape the `n-th` vector by `reshape(n, (ntuple(Returns(1), n-1)..., :))`

2 Likes

Got it ``````julia> outer_product2(vs) =  .*([reshape(vs[d], (ntuple(Returns(1), d-1)..., :)) for d in 1:length(vs)]... )

julia> @benchmark outer_product(\$vs)
BenchmarkTools.Trial: 10000 samples with 5 evaluations.
Range (min … max):  6.045 μs …  1.440 ms  ┊ GC (min … max): 0.00% … 98.98%
Time  (median):     6.857 μs              ┊ GC (median):    0.00%
Time  (mean ± σ):   7.907 μs ± 31.443 μs  ┊ GC (mean ± σ):  8.84% ±  2.21%

▅▃▄▃▁█
███████▇▅▄▃▃▃▂▃▂▃▂▄▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
6.04 μs        Histogram: frequency by time        15.6 μs <

Memory estimate: 5.11 KiB, allocs estimate: 135.

julia> @benchmark outer_product2(\$vs)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
Range (min … max):  1.632 μs … 591.462 μs  ┊ GC (min … max): 0.00% … 99.03%
Time  (median):     1.880 μs               ┊ GC (median):    0.00%
Time  (mean ± σ):   2.032 μs ±   5.930 μs  ┊ GC (mean ± σ):  2.88% ±  0.99%

██▃▁ ▄▂ ▂ ▁▂
▂▄▇████▇████████▆█▆▄▃▄▃▃▃▃▃▃▄▄▂▂▂▂▂▂▂▂▂▂▂▃▄▃▂▂▂▃▂▂▁▁▁▁▁▂▂▁▁ ▃
1.63 μs         Histogram: frequency by time        2.82 μs <

Memory estimate: 752 bytes, allocs estimate: 15.
``````

thanks It’s pretty simple with `Iterators.product`:

``````outer_product3(vs) = Base.splat(*).(Iterators.product(vs...))
``````

Also seems to be faster for this small example:

``````julia> @benchmark outer_product(\$vs)
BenchmarkTools.Trial: 10000 samples with 6 evaluations.
Range (min … max):  5.900 μs … 910.300 μs  ┊ GC (min … max): 0.00% … 98.87%
Time  (median):     6.117 μs               ┊ GC (median):    0.00%
Time  (mean ± σ):   6.828 μs ±  21.970 μs  ┊ GC (mean ± σ):  7.83% ±  2.42%

▆▇█▇▆▆▅▅▄▃▃▂▃▃▂▂▁▂▂▁▂▁ ▁ ▁▁                                ▂
████████████████████████████▇█▇█▆▇▇▆▅▆▅▅▅▅▅▅▄▆▅▅▄▅▃▃▄▂▄▆▅▆▅ █
5.9 μs       Histogram: log(frequency) by time      8.47 μs <

Memory estimate: 5.11 KiB, allocs estimate: 135.

julia> @benchmark outer_product2(\$vs)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
Range (min … max):  1.530 μs … 488.580 μs  ┊ GC (min … max): 0.00% … 98.97%
Time  (median):     1.610 μs               ┊ GC (median):    0.00%
Time  (mean ± σ):   1.756 μs ±   4.889 μs  ┊ GC (mean ± σ):  2.75% ±  0.99%

▅█
▄███▆▃▃▄▅▅▃▃▂▃▄▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▁▁▂▂▂▁▂▂▁▂▂▂▂▂▁▂ ▃
1.53 μs         Histogram: frequency by time        2.96 μs <

Memory estimate: 752 bytes, allocs estimate: 15.

julia> @benchmark outer_product3(\$vs)
BenchmarkTools.Trial: 10000 samples with 176 evaluations.
Range (min … max):  605.114 ns …  13.131 μs  ┊ GC (min … max): 0.00% … 92.48%
Time  (median):     652.841 ns               ┊ GC (median):    0.00%
Time  (mean ± σ):   707.149 ns ± 606.302 ns  ┊ GC (mean ± σ):  4.17% ±  4.69%

█▅▃ ▇▆▅▃▂▂▂▁ ▁   ▂▂▂▁   ▁▂▁▁                                  ▂
█████████████████████▆▇█████▇▇▇▆▆▆▇▆▅▅▆▆▆▁▄▃▄▇█▅▅▅██▅▇▆▆▅▄▆▅▆ █
605 ns        Histogram: log(frequency) by time       1.28 μs <

Memory estimate: 672 bytes, allocs estimate: 7.
``````
2 Likes

for length of vs greater than 3 this could perform better

``````mapreduce(d->reshape(vs[d], (ntuple((x)->1, d-1)..., :)),.*,2:length(vs), init=vs);
``````
1 Like

Nice one.
one can save Base.splat

``````outer_product31(vs) = prod.(Iterators.product(vs...))
``````
3 Likes
``````julia> VS=[rand(rand(2:5)) for _ in 1:10]
10-element Vector{Vector{Float64}}:
[0.7809510482067282, 0.9707013914054325, 0.3620120022634721, 0.7488162372721611, 0.2208620150439493]
[0.6110369641230907, 0.5098682950578028, 0.5049683337707953, 0.8757125530915953]
[0.9385327815487692, 0.46337103561713033, 0.7475564936185806]
[0.12816029157324627, 0.9877131460915131, 0.7182428947165578]
[0.6111863554757773, 0.2318165661824203, 0.44984729314945215]
[0.9945688985718182, 0.9665853587265798, 0.6521171926682267, 0.29032509439145715]
[0.60695142806559, 0.2895700161283634]
[0.7383875082243384, 0.25652754729435023, 0.7715442370365116, 0.6833019416092574, 0.21299555398838987]
[0.7822395058310352, 0.520211099994906, 0.8233746349481897, 0.26566768583948375]
[0.37286173425537894, 0.8832868074013911, 0.18649640224567243, 0.3432979876570885]

julia> outer_product2(vs) =  .*([reshape(vs[d], (ntuple(Returns(1), d-1)..., :)) for d in 1:length(vs)]... )
outer_product2 (generic function with 1 method)

julia> outer_product4(vs)=mapreduce(d->reshape(vs[d], (ntuple((x)->1, d-1)..., :)),.*,2:length(vs), init=vs);

julia> outer_product3(vs) = Base.splat(*).(Iterators.product(vs...))
outer_product3 (generic function with 1 method)

julia> outer_product31(vs) = prod.(Iterators.product(vs...))
outer_product31 (generic function with 1 method)

julia> outer_product5(vs)=reduce((l,r)->SplitApplyCombine.product(*, l,r), VS);

julia> @btime outer_product2(VS);
3.562 ms (63 allocations: 2.64 MiB)

julia> @btime outer_product3(VS);
7.056 ms (10 allocations: 29.00 MiB)

julia> @btime outer_product31(VS);
7.170 ms (10 allocations: 29.00 MiB)

julia> @btime outer_product4(VS);
1.459 ms (50 allocations: 3.52 MiB)

julia> @btime outer_product5(VS);
1.383 ms (129 allocations: 3.53 MiB)
``````

Note that broadcasting over an iterator will collect it first; `map` avoids this:

``````julia> @btime outer_product31(\$VS);  # prod.(Iterators.product(vs...))
min 3.102 ms, mean 4.220 ms (10 allocations, 19.34 MiB)

julia> @btime outer_product4(\$VS);
min 1.545 ms, mean 1.882 ms (100 allocations, 2.40 MiB)

julia> outer_product32(vs) = map(prod, Iterators.product(vs...));

julia> @btime outer_product32(\$VS);
min 296.000 μs, mean 463.584 μs (7 allocations, 1.76 MiB)

julia> VS |> summary
"10-element Vector{Vector{Float64}}"
``````

None of these are type-stable, as `ndims` of the final array depends on a length. Dealing with a tuple instead of a vector may be much better.

``````julia> @btime outer_product32(\$vs);  # splat within the timing
min 807.411 ns, mean 847.781 ns (5 allocations, 288 bytes)

julia> @btime outer_product32(\$(Tuple(vs)));  # make Tuple before timing
min 83.289 ns, mean 97.792 ns (1 allocation, 160 bytes)

julia> vs |> summary  # smaller inputs than VS
"3-element Vector{Vector{Float64}}"
``````
3 Likes

Maybe then this topic should be the start of a dedicated function in Julia  (where I think that the inputs should be the various vectors rather than a vector of vector)

`SplitApplyCombine.product` and `productview` take two containers and do the outer product

``````julia> product(+, [1,2], [1,2,3])
2×3 Array{Int64,2}:
2  3  4
3  4  5

julia> productview(+, [1,2], [1,2,3])
2×3 ProductArray{Int64,2,typeof(+),Array{Int64,1},Array{Int64,1}}:
2  3  4
3  4  5
``````

It could be useful to make these take more than two container arguments.

``````Base.splat(+).(Base.product([1,2], [1,2,3], [-1, -2]))
``````
``````using SplitApplyCombine
foldl( (x,y)->product(+,x,y), [[1,2],[1,2,3],[-1,-2]])
``````
1 Like

One could also reshape a `kron`-product:

``````outer(vs) = reshape(kron(reverse(vs)...), length.(vs)...)
``````

This is the same as your `outer_product`.

Agree!

``````outer(vs...) = reshape(kron(reverse(vs)...), length.(vs))
``````

Julia gets dinged around minute 15 since broadcasting does not equate directly with an outer product in all circumstances. This gets explained further around minute 32 where he takes issues with `transpose` on `Vector`. See also this prior post:

ThIs thread clearly shows that a true outer product is also accessible in Julia.

Hoekstra’s original version without using `I` is as follows:

``````# check fo see if matrix is an X matrix
# where only the main diagonals are non-zero
function checkmatrix(grid)
n = size(grid, 1)
i = 1:n .== 1:n |> transpose
min.(grid, 1) == max.(i, reverse(i, dims=1))
end
``````

Using `Iterators.product`, `map`, and `allequal` we could do the following.

``````function checkmatrix2(grid)
n = size(grid, 1)
i = map(allequal, Iterators.product(1:n, 1:n))
min.(grid, 1) == max.(i, reverse(i, dims=1))
end
``````

nice exercise.

some alternatives, just to ask a few questions.

Why, using the following expression, does product not get a Matrix{Bool}?

``````using SplitApplyCombine
julia> i=product((x,y)->x==y || x+y==(n+1),1:n,1:n)
9×9 Matrix{Any}:
true  false  false  false  false  false  false  false   true
false   true  false  false  false  false  false   true  false
false  false   true  false  false  false   true  false  false
false  false  false   true  false   true  false  false  false
false  false  false  false   true  false  false  false  false
false  false  false   true  false   true  false  false  false
false  false   true  false  false  false   true  false  false
false   true  false  false  false  false  false   true  false
true  false  false  false  false  false  false  false   true

julia> i=identity.(product((x,y)->x==y || x+y==(n+1),1:n,1:n))
9×9 BitMatrix:
1  0  0  0  0  0  0  0  1
0  1  0  0  0  0  0  1  0
0  0  1  0  0  0  1  0  0
0  0  0  1  0  1  0  0  0
0  0  0  0  1  0  0  0  0
0  0  0  1  0  1  0  0  0
0  0  1  0  0  0  1  0  0
0  1  0  0  0  0  0  1  0
1  0  0  0  0  0  0  0  1
``````
``````function checkmatrix3(grid)
n = size(grid, 1)
findall(map(!=(0.), grid).!= product((x,y)->x==y || x+y==(n+1),1:n,1:n))
end
``````

Now the main question. Probably my doubt comes from the poor knowledge of English, but I wonder if a matrix of all zeros does not logically satisfy the condition.

The expression “a matrix where only the main diagonals are non-zero” implies that the non-diagonally principal elements are all zero, but it says nothing about the diagonally principal elements. Is that it?