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
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}}"
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)
I was thinking about this after viewing the following video by Conor Hoekstra:
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
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?
