# Outer product of more than two vectors

**URL:** <https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310>\
**Category:** General Usage\
**Tags:** question, performance, linearalgebra\
**Created:** [March 19, 2023, 3:48pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310 "2023-03-19T15:48:10Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 19, 2023, 3:48pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/1 "2023-03-19T15:48:11Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [March 19, 2023, 3:52pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/2 "2023-03-19T15:52:23Z")

</div>

```julia
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

```

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 19, 2023, 3:53pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/3 "2023-03-19T15:53:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [March 19, 2023, 3:55pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/4 "2023-03-19T15:55:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 19, 2023, 4:13pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/5 "2023-03-19T16:13:48Z")

</div>

Got it 🙂

```julia
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 🙂

---

<div class="post-metadata">

**Author:** ![yha](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yha/32/3502_2.png) [@yha](https://discourse.julialang.org/u/yha)\
**Post date:** [March 19, 2023, 4:30pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/6 "2023-03-19T16:30:20Z")

</div>

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

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

```

Also seems to be faster for this small example:

```julia
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.

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [March 19, 2023, 5:36pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/7 "2023-03-19T17:36:35Z")

</div>

for length of vs greater than 3 this could perform better

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

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [March 19, 2023, 5:43pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/8 "2023-03-19T17:43:49Z")

</div>

Nice one.  
one can save Base.splat

```julia
outer_product31(vs) = prod.(Iterators.product(vs...))

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [March 19, 2023, 6:17pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/9 "2023-03-19T18:17:54Z")

</div>

```julia
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[1]);

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)

```

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [March 19, 2023, 6:38pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/10 "2023-03-19T18:38:55Z")

</div>

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

```julia
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
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}}"

```

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [March 20, 2023, 8:10am UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/12 "2023-03-20T08:10:18Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [March 25, 2023, 10:30pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/13 "2023-03-25T22:30:47Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [March 25, 2023, 10:55pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/14 "2023-03-25T22:55:16Z")

</div>

```julia
Base.splat(+).(Base.product([1,2], [1,2,3], [-1, -2]))

```

```julia
using SplitApplyCombine
foldl( (x,y)->product(+,x,y), [[1,2],[1,2,3],[-1,-2]])

```

---

<div class="post-metadata">

**Author:** ![Jollywatt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jollywatt/32/202198_2.png) [@Jollywatt](https://discourse.julialang.org/u/Jollywatt)\
**Post date:** [March 25, 2023, 11:41pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/15 "2023-03-25T23:41:53Z")

</div>

> [@sylvaticus](#):
>
> Is there a built-in function for an outer product of more than two vectors

One could also reshape a `kron`-product:

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

```

This is the same as your `outer_product`.

---

<div class="post-metadata">

**Author:** ![Jollywatt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jollywatt/32/202198_2.png) [@Jollywatt](https://discourse.julialang.org/u/Jollywatt)\
**Post date:** [March 25, 2023, 11:48pm UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/16 "2023-03-25T23:48:37Z")

</div>

> [@sylvaticus](#):
>
> (where I think that the inputs should be the various vectors rather than a vector of vector)

Agree!

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

```

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [March 29, 2023, 6:31am UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/17 "2023-03-29T06:31:19Z")

</div>

I was thinking about this after viewing the following video by Conor Hoekstra:

[![](https://global.discourse-cdn.com/julialang/original/3X/2/b/2b46fc4cd37dff595e491e9a8a642d23ce24babd.jpeg "APL vs BQN vs J vs Q vs NumPy vs Julia vs R") ](https://www.youtube.com/watch?v=8ynsN4nJxzU)

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:

> [@Should transpose of a Vector be a Matrix?](https://discourse.julialang.org/t/should-transpose-of-a-vector-be-a-matrix/87444):
>
> As Julia has a true 1D vector data type, shouldn’t the transpose of a Vector just be a Vector (essentially an identity operation)? See Conor Hoekstra’s comments in [APL vs BQN vs J vs Q vs NumPy vs Julia vs R - YouTube](https://youtu.be/8ynsN4nJxzU?t=1947) I assume this stems from MATLAB where a vector is actually a 2D array and that you can, of course, transpose. What are your thoughts about this?

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

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

```julia
# 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.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [March 29, 2023, 7:57am UTC](https://discourse.julialang.org/t/outer-product-of-more-than-two-vectors/96310/18 "2023-03-29T07:57:17Z")

</div>

nice exercise.

some alternatives, just to ask a few questions.

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

```julia
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

```

```julia
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?