Cross Product in Julia

Viewed 173

Consider the Matrix

A = [1., 2., 3.; 4., 5., 6]

and the Vector

v = [0.1, 0.2, 0.4]

What is the most efficient "Julianic" way to obtain a new Matrix whose rows are [cross(A[1, :], v); cross([2, :], v)] ? The result for the numbers above should be [0.2, -0.1, 0.; 0.8, -0.1, 0.3]

How to generalize for N rows in A?

4 Answers

Try using the cross-product operator matrix

cross_op(v) = [0.0 -v[3] v[2]; v[3] 0.0 -v[1]; -v[2] v[1] 0.0];

and then use it as

w = A*cross_op(v)

You can do it with dot broadcasting like this:

A = [1 2 3; 4 5 6];

v = [0.1, 0.2, 0.4];

cross.(eachrow(A), Ref(v))
2-element Vector{Vector{Float64}}:
 [0.19999999999999996, -0.09999999999999998, 0.0]
 [0.7999999999999998, -1.0, 0.30000000000000004]

Note that the original A in the question is not valid syntax. In Julia, a literal matrix is written with rows having spaces between elements.

The obvious built-in function for this is mapslices, not mentioned so far. But unfortunately it is slow, until Julia 1.9 comes out.

Julia 1.9 will also let you replace most uses of the fast (but hard-to-discover) reduce(hcat, ...) with stack. You can use this today via the Compat.jl package, which contains features from the future.

julia> A = [1. 2. 3.; 4. 5. 6]  # note that the original had extra commas
2×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0

julia> v = [.1, .2, .4];

julia> B = @btime mapslices(x -> cross(x, $v), $A; dims=2)
  min 2.866 μs, mean 3.137 μs (44 allocations, 2.08 KiB)      # Julia 1.8
  min 244.659 ns, mean 305.486 ns (4 allocations, 352 bytes)  # Julia 1.9
2×3 Matrix{Float64}:
 0.2  -0.1  0.0
 0.8  -1.0  0.3

julia> using Compat  # loads some forthcoming features

julia> B == @btime stack(cross(x, $v) for x in eachrow($A); dims=1)
  min 129.778 ns, mean 141.398 ns (3 allocations, 272 bytes)
true

julia> B == @btime transpose(reduce(hcat, [cross(a, $v) for a in eachrow($A)]))
  min 169.444 ns, mean 190.345 ns (4 allocations, 336 bytes)
true  # note this returns a lazy Transpose, not a Matrix

As others have said, using columns not rows is probably more natural, if you can:

julia> A2 = [1. 4.; 2. 5.; 3. 6.];

julia> B2 = mapslices(x -> cross(x, v), A2; dims=1)

julia> B2 == @btime stack(cross(x, $v) for x in eachcol($A2));
  min 116.812 ns, mean 129.633 ns (3 allocations, 272 bytes)

julia> B2 == @btime reduce(hcat, [cross(a, $v) for a in eachcol($A2)])
  min 168.212 ns, mean 189.909 ns (4 allocations, 336 bytes)
true

For handling many short vectors, you may also want to investigate using StaticArrays.jl:

julia> using StaticArrays

julia> v3 = SA[.1, .2, .4];

julia> A3 = reinterpret(reshape, typeof(v3), A2)  # view each col as a vector
2-element reinterpret(reshape, SVector{3, Float64}, ::Matrix{Float64}) with eltype SVector{3, Float64}:
 [1.0, 2.0, 3.0]
 [4.0, 5.0, 6.0]

julia> B2 == @btime stack(cross(x, $v3) for x in $A3);
  min 61.735 ns, mean 65.846 ns (1 allocation, 112 bytes)

When expecting the result to be the same type as A create it using similar and fill it in a loop with the result of cross.

A = [1. 2. 3.; 4. 5. 6]
v = [.1, .2, .4]

using LinearAlgebra  #for command cross

B = similar(A)
for (i, x) in enumerate(eachrow(A))
    B[i, :] = cross(x, v)
end
B
#2×3 Matrix{Float64}:
# 0.2  -0.1  0.0
# 0.8  -1.0  0.3

using BenchmarkTools
@btime begin
    B = similar($A)
    for (i, x) in enumerate(eachrow($A))
        B[i, :] = cross(x, $v)
    end
end
#  78.340 ns (3 allocations: 272 bytes)

As pointed out by @DNF it will be faster to go along columns instead of rows.

A2 = [1. 4.; 2. 5.; 3. 6.] #Changed orientation of A

@btime begin
    B = similar($A2)
    for (j, x) in enumerate(eachcol($A2))
        B[:, j] = cross(x, $v)
    end
end
#   77.214 ns (3 allocations: 272 bytes)

Another option will be to use hcat in reduce.

@btime reduce(hcat, [cross(a, $v) for a in eachrow($A)])
#  103.328 ns (4 allocations: 336 bytes)

Compared with the similar method from @DNF

@btime hcat([cross(a, $v) for a in eachrow($A)]...)
#  124.479 ns (4 allocations: 336 bytes)

And to those from other answers

@btime cross.(eachrow($A), Ref($v))  #@AboAmmar - but returns a vector
#  104.453 ns (5 allocations: 384 bytes)

@btime vcat([permutedims(cross(r,$v)) for r in eachrow($A)]...)  #@Dan Getz
#  180.323 ns (8 allocations: 528 bytes)

@btime $A * [0.0 -$v[3] $v[2]; $v[3] 0.0 -$v[1]; -$v[2] $v[1] 0.0]  #@John Alexiou
#  263.570 ns (9 allocations: 416 bytes)

@btime [cross($A[1, :], $v), cross($A[2, :], $v)]  #@JustLearning - but returns a vector
#  107.073 ns (5 allocations: 384 bytes)

@btime mapslices(x -> cross(x, $v), $A; dims=2)  #@mcabbott
#  1.887 μs (44 allocations: 2.08 KiB)

@btime [cross(a, $v) for a in eachrow($A)]  #but returns a vector
#  73.807 ns (3 allocations: 224 bytes)

@btime [cross(a, $v) for a in eachcol($A2)]  #but returns a vector
#  71.950 ns (3 allocations: 224 bytes)

@mcabbott pointed that in 1.9 stack will be available which is also fast, especially when A can be used by col.

using Compat #for stack from v1.9
@btime stack(cross(x, $v) for x in eachrow($A); dims=1)
#  93.161 ns (3 allocations: 272 bytes)

@btime stack(cross(x, $v) for x in eachcol($A2))
#  79.622 ns (3 allocations: 272 bytes)

@mcabbott also showed that when using StaticArrays it could be done even faster.

using StaticArrays
v3 = SA[.1, .2, .4];
A3 = reinterpret(reshape, typeof(v3), A2)
@btime stack(cross(x, $v3) for x in $A3)
#  43.051 ns (1 allocation: 112 bytes)

What could also be done with a loop.

@btime begin
    B = similar($A2)
    for (i, x) in enumerate($A3)
        B[:, i] = cross(x, $v3)
    end
end
#  40.591 ns (1 allocation: 112 bytes)

@btime [cross(x, $v3) for x in $A3]  #But returns a vector
#  36.914 ns (1 allocation: 112 bytes)

Currently it look like that the most efficient "Julianic" way to obtain a new Matrix is a for loop in case the given data should be used as it is. In case it could be stored in a different way stack is as efficient as the loop but more convenient. Using StaticArrays will improve the performance of both.
But if the matrix has many rows

A = randn((10000, 3))

@btime begin
    B = similar($A)
    for (i, x) in enumerate(eachrow($A))
        B[i, :] = cross(x, $v)
    end
end
#  246.372 μs (10002 allocations: 1015.67 KiB)

@btime $A * [0.0 -$v[3] $v[2]; $v[3] 0.0 -$v[1]; -$v[2] $v[1] 0.0]
#  30.958 μs (10 allocations: 234.72 KiB)

using StaticArrays
v3 = SA[.1, .2, .4];
A3 = reinterpret(reshape, typeof(v3), A')
@btime stack(cross(x, $v3) for x in $A3)
  96.000 μs (2 allocations: 234.42 KiB)

the cross-product operator matrix given in the answer of @John Alexiou is the fastest!

Related