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

**URL:** <https://discourse.julialang.org/t/vectors-matrices-with-units-vs-units-for-individual-entries-unitful-jl/50453>\
**Category:** General Usage\
**Tags:** unitful\
**Created:** [November 19, 2020, 6:31pm UTC](https://discourse.julialang.org/t/vectors-matrices-with-units-vs-units-for-individual-entries-unitful-jl/50453 "2020-11-19T18:31:15Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![airpmb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/airpmb/32/7826_2.png) [@airpmb](https://discourse.julialang.org/u/airpmb)\
**Post date:** [November 19, 2020, 6:31pm UTC](https://discourse.julialang.org/t/vectors-matrices-with-units-vs-units-for-individual-entries-unitful-jl/50453/1 "2020-11-19T18:31:15Z")

</div>

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):

```julia
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
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
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:

```julia
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
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
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.

---

<div class="post-metadata">

**Author:** ![FedericoStra](https://avatars.discourse-cdn.com/v4/letter/f/76d3ee/32.png) [@FedericoStra](https://discourse.julialang.org/u/FedericoStra)\
**Post date:** [November 20, 2020, 8:51pm UTC](https://discourse.julialang.org/t/vectors-matrices-with-units-vs-units-for-individual-entries-unitful-jl/50453/2 "2020-11-20T20:51:38Z")

</div>

> [@airpmb](#):
>
> ```julia
> 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
> 
> ```

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`:

```julia
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

```

* * *

> [@airpmb](#):
>
> surprisingly (or maybe not) specifying units on each entry rather than for the whole array seems significantly faster:

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](http://painterqubits.github.io/Unitful.jl/stable/highlights/#Arrays-1)).

```julia
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.

---

<div class="post-metadata">

**Author:** ![jlawrie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlawrie/32/26283_2.png) [@jlawrie](https://discourse.julialang.org/u/jlawrie)\
**Post date:** [August 19, 2022, 11:35am UTC](https://discourse.julialang.org/t/vectors-matrices-with-units-vs-units-for-individual-entries-unitful-jl/50453/3 "2022-08-19T11:35:42Z")

</div>

Hello -

thanks for the notes, I found these helpful whilst troubleshooting a different possible bug with Unitful.jl!

Just for completeness’ sake, I noticed that another way (see `f_vec1()` below) of constructing a vector of “unitful” values is even ‘worse’ for performance (and at the moment it’s this implementation which is causing my other problem, unfortunately).

```julia
julia> f_vec1() = Quantity{Number}[i*u"L/minute" for i in 0:100]
f_vec1 (generic function with 1 method)

julia> f_vec2() = [i*u"L/minute" for i in 0:100]
f_vec2 (generic function with 1 method)

julia> f_vec3() = [i for i in 0:100]u"L/minute"
f_vec3 (generic function with 1 method)

julia> @btime f_vec1();

  1.196 μs (102 allocations: 2.45 KiB)

julia> @btime f_vec2();

  109.337 ns (1 allocation: 896 bytes)

julia> @btime f_vec3();

  250.504 ns (2 allocations: 1.75 KiB)

```

Thanks!

* * *

Julia Version:

```julia
julia> versioninfo()
Julia Version 1.8.0
Commit 5544a0fab76 (2022-08-17 13:38 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin21.4.0)
  CPU: 4 × Intel(R) Core(TM) i5-5250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, broadwell)
  Threads: 1 on 2 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 

```

Unitful.jl Version:

```julia
(@v1.8) pkg> st Unitful
Status `~/.julia/environments/v1.8/Project.toml`
  [1986cc42] Unitful v1.11.0

```

BenchmarkTools.jl Version:

```julia
(@v1.8) pkg> st BenchmarkTools
Status `~/.julia/environments/v1.8/Project.toml`
  [6e4b80f9] BenchmarkTools v1.3.1

```

* * *