Type instability in array concatenation

General Usage Performance
1

I have some simple code to make a left-quaternion multiplication matrix:

function L_q(q::Vector{T}) where {T<:Real}
    return [q[1] -q[2:4]'; q[2:4] q[1] * I + SkewSymmetricMatrix(q[2:4])]
end

But the problem is that the hvcat operation makes the return type unstable. Here’s the output from running

function L_q(q::Vector{T}) where {T<:Real}
    return [q[1] -q[2:4]'; q[2:4] q[1] * I + SkewSymmetricMatrix(q[2:4])]
end
function SkewSymmetricMatrix(q::Vector{T}) where {T<:Real}
    return [T(0) -q[3] q[2]; q[3] T(0) -q[1]; -q[2] q[1] T(0)]
end

q = [1.0, 0.0, 0.0, 0.0]
@code_warntype L_q(q)

The output:

Variables
  #self#::Core.Compiler.Const(L_q, false)
  q::Array{Float64,1}

Body::Any
1 ─ %1  = Core.tuple(2, 2)::Core.Compiler.Const((2, 2), false)
│   %2  = Base.getindex(q, 1)::Float64
│   %3  = (2:4)::Core.Compiler.Const(2:4, false)
│   %4  = Base.getindex(q, %3)::Array{Float64,1}
│   %5  = Base.adjoint(%4)::Adjoint{Float64,Array{Float64,1}}
│   %6  = -%5::Adjoint{Float64,Array{Float64,1}}
│   %7  = (2:4)::Core.Compiler.Const(2:4, false)
│   %8  = Base.getindex(q, %7)::Array{Float64,1}
│   %9  = Base.getindex(q, 1)::Float64
│   %10 = (%9 * Main.I)::UniformScaling{Float64}
│   %11 = (2:4)::Core.Compiler.Const(2:4, false)
│   %12 = Base.getindex(q, %11)::Array{Float64,1}
│   %13 = Main.SkewSymmetricMatrix(%12)::Array{Float64,2}
│   %14 = (%10 + %13)::Array{Float64,2}
│   %15 = Base.hvcat(%1, %2, %6, %8, %14)::Any
└──       return %15

There’s an ::Any on line %15. Is the problem because it’s concatenating an adjoint with other arrays? I read that I could use copy(…) but just doing copy(q[2:4]') didn’t work.

Would appreciate any ideas to fix this. Thanks!

2

Please provide an MWE.

3

Last time I looked at the implementation of hvcat (mind you, it was years ago, so things may have changed) it was wonderfully general but inherently type unstable. A better chance for type stability would be to use a typed hvcat,

function L_q(q::Vector{T}) where {T<:Real}
    return T[q[1] -q[2:4]'; q[2:4] q[1] * I + SkewSymmetricMatrix(q[2:4])]
end

but it doesn’t seem to be enough. Adding a return type declaration does help with the return type though:

function L_q(q::Vector{T})::Matrix{T} where {T<:Real}
    return T[q[1] -q[2:4]'; q[2:4] q[1] * I + SkewSymmetricMatrix(q[2:4])]
end
5

Thanks for the suggestion! I wonder why typed_hvcat doesn’t fix the issue. Is there documentation for this function anywhere?