Output type of binary operations involving StaticArrays

General Usage Performance
1

binary operations on a StaticArray and Array give outputs with quite-inconsistent types. For example, intuitively SVector .+ Vector should give a SVector rather than a Vector?

julia> using StaticArrays

julia> typeof(SVector{3, Float64}(randn(3)) .+ randn(3) )
Array{Float64,1}

julia> typeof(SVector{3, Float64}(randn(3)) .* randn(3) )
Array{Float64,1}

julia> typeof(SVector{3, Float64}(randn(3)) .+ 1.0 )
SArray{Tuple{3},Float64,1,3}

julia> typeof(SVector{3, Float64}(randn(3)) .* 2.0 )
SArray{Tuple{3},Float64,1,3}

julia> typeof(SMatrix{3, 2, Float64}(randn(3, 2)) * randn(2) )
SArray{Tuple{3},Float64,1,3}

julia> typeof(randn(1, 3) * SMatrix{3, 2, Float64}(randn(3, 2)) )
Array{Float64,2}

how does the output type be determined? could I override them? thanks.

2

I am not sure about this — both options can be reasonable and consistent.

If you want static results, use statically sized decorators (see the docs):

using StaticArrays
a = SVector{3, Float64}(randn(3))
b = SizedVector{3}(randn(3))
a .+ b # 3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3)
3

thanks.
what I mean by inconsistency is some op.(::StaticArray, ::Array) give a StaticArray while other op(). gives Array.
that means we have to test every op(). before knowing the type it gives.

btw, does SizedVector{N}() cause any runtime?

actually, how could I control the output type of generally op.(::A, ::B)? is it about broadcasting? thanks.

4

I am not sure I understand why you think this. It may be better to just write generic code, and control your inputs when it matters.

Quite the opposite, the purpose is to aid inference. Again, please see the StaticArrays manual.

5

seems like I had some confusions about x .+ y vs. x .+ 1.0. To clarify, I try the following:

module Testing

import Base.size, Base.getindex
export A, B

struct A{T <: AbstractVector{Float64} } <: AbstractVector{Float64}
    data::T
end

struct B{T <: AbstractVector{Float64} } <: AbstractVector{Float64}
    data::T
end

Base.size(x::A) = size(x.data)
Base.size(x::B) = size(x.data)

Base.getindex(x::A, i) = x.data[i]
Base.getindex(x::B, i) = x.data[i]

end

using Main.Testing
a = A(randn(3) )
b = B(randn(3) )

julia> typeof(exp.(a) )
Array{Float64,1}
julia> typeof(a .+ a)
Array{Float64,1}
julia> typeof(a .+ b)
Array{Float64,1}
julia> typeof(b .+ a)
Array{Float64,1}
julia> typeof(a .+ 1.0)
Array{Float64,1}
julia> typeof(b .+ 1.0)
Array{Float64,1}

the last five typeof() results show that, by default, exp.(::A), .+(::A, ::B) and .+(::A, ::Float64) always output a Vector (but not a A nor a B).

So now I know that .+(::SVector, ::Vector) giving a Vector is just the default behavior.

However, there’re important exceptions e.g. exp.(::SVector), .+(::SVector, ::Float64) and .+(::SVector, ::SizedVector): they output SVector but not the default Vector !

How could it be done :dizzy_face:? An example making exp.(::A), .+(::A, ::B) and .+(::A, ::Float32) outputting a A is highly appreciated, thanks.

6

I would agree with that… in an ideal world. The broadcasting machinery is sufficiently intricate that it might not be possible/reasonable though. File an issue?

7

In some cases wrapping in SizedVector may slow things down since one additional allocation happens (at least in stable versions of Julia).

This should be possible but at the cost of making custom broadcasting in StaticArrays more complicated. It’d have to figure out which arrays should have sizes statically checked and which should have run-time checks.

1 Like
8

That’s interesting to know, can you please give an example?

9

For example:

julia> f(n) = SizedVector{2}([n, n+1])
f (generic function with 1 method)

julia> g(n) = [n, n+1]
g (generic function with 1 method)

julia> @btime f(10)
  37.265 ns (2 allocations: 112 bytes)
2-element SizedArray{Tuple{2},Int64,1,1} with indices SOneTo(2):
 10
 11

julia> @btime g(10)
  28.605 ns (1 allocation: 96 bytes)
2-element Array{Int64,1}:
 10
 11

Although the worse performance is limited to very specific circumstances (lots of allocations that can’t be optimized away). similar for SizedArray was also surprisingly slow until a few months ago: Faster similar for SizedArray by mateuszbaran · Pull Request #625 · JuliaArrays/StaticArrays.jl · GitHub

1 Like
10

Thanks. My suggestion was more about the original example of @tomtom, something like

julia> using StaticArrays, BenchmarkTools

julia> f(x, n) = x .+ SizedVector{2}([n, n+1])
f (generic function with 1 method)

julia> g(x, n) = x .+ [n, n+1]
g (generic function with 1 method)

julia> x = SVector(3, 4)
2-element SArray{Tuple{2},Int64,1,2} with indices SOneTo(2):
 3
 4

julia> @btime f($x, $10)
  33.359 ns (1 allocation: 96 bytes)
2-element SArray{Tuple{2},Int64,1,2} with indices SOneTo(2):
 13
 15

julia> @btime g($x, $10)
  62.608 ns (2 allocations: 192 bytes)
2-element Array{Int64,1}:
 13
 15
1 Like
11

more examples of “inconsistency”:

julia> exp.(SizedVector{3, Float64}(randn(3) ) )
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):

julia> 1.0 .+ SizedVector{3, Float64}(randn(3) )
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):

operations on SizedVector gives SVector rather than SizedVector !

but now SizedVector “loses” to SVector or Vector:

julia> SVector{3, Float64}(randn(3) ) .+ SizedVector{3, Float64}(randn(3) )
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):

julia> (1.0:3.0) .+ SizedVector{3, Float64}(randn(3) )
3-element Array{Float64,1}:

that said, it’s hard to predict the type of output unless a testing is done.

12

As far as I understand StaticArrays the general rule is that when all arguments to a broadcast are StaticArrays or scalars (or some wrapped StaticArrays), then SArray will be returned (because it’s the fastest thing). You’ll see it in many other places as well, StaticArrays operations usually return SArray when they can.

1 Like