Vectors/matrices with units vs. units for individual entries (Unitful.jl)

I ran into something interesting yesterday - I wrote this function to create the cross-product matrix for 3-vector (sometimes known as the tilde operator):

tilde(x) = [0 -x[3] x[2]; x[3] 0 -x[1]; -x[2] x[1] 0]

If I then use it with an input that has units, it seemingly works OK:

julia> v = [0.1 0.2 0.3]u"rad/s"
1×3 Array{Quantity{Float64,𝐓⁻¹,Unitful.FreeUnits{(rad, s⁻¹),𝐓⁻¹,nothing}},2}:
 0.1 rad s⁻¹  0.2 rad s⁻¹  0.3 rad s⁻¹
julia> ṽ = tilde(v)
3×3 Array{Quantity{Float64,D,U} where U where D,2}:
          0.0  -0.3 rad s⁻¹   0.2 rad s⁻¹
  0.3 rad s⁻¹           0.0  -0.1 rad s⁻¹
 -0.2 rad s⁻¹   0.1 rad s⁻¹           0.0

Trouble begins, though, when I try to use this matrix:

julia> ṽ*ṽ'
ERROR: DimensionError: 0.0 and 0.09 rad² s⁻² are not dimensionally compatible.
Stacktrace:
 [1] +(::Quantity{Float64,NoDims,Unitful.FreeUnits{(),NoDims,nothing}}, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at /Users/patrick/.julia/packages/Unitful/1t88N/src/quantities.jl:137
 [2] +(::Float64, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at ./promotion.jl:311
 [3] +(::Float64, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at ./operators.jl:538
 [4] matmul3x3!(::Array{Any,2}, ::Char, ::Char, ::Array{Quantity{Float64,D,U} where U where D,2}, ::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:996
 [5] generic_matmatmul!(::Array{Any,2}, ::Char, ::Char, ::Array{Quantity{Float64,D,U} where U where D,2}, ::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:726
 [6] mul! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:411 [inlined]
 [7] mul! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:208 [inlined]
 [8] *(::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.Adjoint{Any,Array{Quantity{Float64,D,U} where U where D,2}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:153
 [9] top-level scope at REPL[10]:100:

I didn’t want to write the function in a way that requires Unitful.jl and found that this does the trick:

function tilde(x)
    z = zero(x[1]) # trick to make this play nice with Unitful.jl
    [z -x[3] x[2]; x[3] z -x[1]; -x[2] x[1] z]
end

The output of this has the same units for every entry and so multiplies without error.

I guess the general question this raises to me is whether there’s a best practice in general for writing a function like mine so that it will play nice with Unitful (and possibly other packages that extend quantities in Julia in interesting ways).

This also got me thinking about whether there’s a way for a vector or matrix as a whole to have units rather than each entry individually having its own units. This does indeed seem to be the case as the type of the array is different when every entry has the same units (or if the units are specified for the array as a whole):

julia> [1u"m" 1u"s"] # entries with different units
1×2 Array{Quantity{Int64,D,U} where U where D,2}:
 1 m  1 s

julia> [1u"m" 1u"m"] # all entries same units
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
 1 m  1 m

julia> [1 1]u"m" # entries unit-less, array given units
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
 1 m  1 m

I’m not sure this can help in my example given the way the matrix is constructed – In the function, how could Julia know that what I want is a matrix with all entries the same type as those of x? Only, I guess, if the code that constructs such a matrix were to assume any zero entry has the same units as the nonzero entries, whenever the latter have the same units for each entry.

Incidentally I wondered if the two ways of creating the “homogeneous units” array were the same under the hood; surprisingly (or maybe not) specifying units on each entry rather than for the whole array seems significantly faster:

julia> fa() = [1u"m" 2u"m"]
fa (generic function with 1 method)

julia> fb() = [1 2]u"m"
fb (generic function with 1 method)

julia> fa() == fb()
true

julia> @btime fa()
  29.271 ns (1 allocation: 96 bytes)
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
 1 m  2 m

julia> @btime fb()
  69.644 ns (2 allocations: 192 bytes)
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
 1 m  2 m

Also relatedly: Obviously for 3x3 cross product matrices there’s not much performance gain to be had but I’m wondering why there isn’t a skew symmetric special matrix type in the LinearAlgebra package.

I think in general it is even better to do z = zero(eltype(x)). Imagine a case where x might be empty, or maybe x[1] does not exist because x is an OffsetArray.

Moreover, in your case you can probably use StaticArrays:

using StaticArrays
function tilde(x)
    z = zero(eltype(x))
    SA[z -x[3] x[2]; x[3] z -x[1]; -x[2] x[1] z]
end

You can debug what the functions fa and fb do with @code_typed. fb does not return a (unitless) array with an associated unit; rather, it construct a unitless array and then multiplies it times a unit, hence producing the same type of object as fa. However, the way this array is constructed is different. You can see with @code_typed that at the end there is a broadcast to perform the multiplication. To get a better sense of the performance, you can check @code_llvm and @code_native too and verify that fb generates substantially more code than fa.


As far as I can tell, there is no way to create arrays with units. The reason is that all the 5 constructors of Quantity require that the quantity in question is a subtype of Number. Still, if the array is homogeneous (every entry has the same unit), then there should be no performance penalty (see also here).

v = [0.1, 0.2, 0.3]u"rad/s"
ṽ = tilde(v)
foo(w) = w * w'
@code_typed foo(ṽ)
@code_llvm foo(ṽ)
@code_native foo(ṽ)

All types are correctly inferred, so the units are effectively tracked at compile time and don’t show up in the code.

1 Like