UnitfulTensors.jl

A package for working efficiently with arrays of physical quantities.

Currently it is experimental and mostly focused on linear algebra.

Features

• Support for arrays with mixed units
• Indexing and assignment: scalar indices, multidimensional indices, logical indexing, omitted/trailing singleton dimensions
• Arithmetic: +, -, *, \, /, ^, sqrt, dot, cross, kron, inv, pinv, lmul!, rmul!, ldiv!, rdiv!
• Transcendental functions of matrices
• adjoint, transpose
• det, logdet, logabsdet, tr
• Factorizations: lu, bunchkaufman, cholesky, ldlt, eigen, hessenberg, schur, svd, lq, qr
• Norm and related functions: cond, condskeel, norm, normalize, nullspace, opnorm, rank
• Rotations and reflections: rotate!, reflect!, givens
• Solving Lyapunov and Sylvester equations
• Support for sparse matrices
• Moderate runtime overhead, often negligible
• Zero-overhead mode with units_off: check the units on a small-scale problem > turn them off > proceed with a large-scale problem

Not implemented yet:

• General array manipulation: hcat, vcat, repeat, reshape, permutedims, sort, etc.
• Elementwise operations

Future plans:

• Add support for frames of reference
• Further optimizations: NoAxisDimensions, InverseAxisDimensions, memoization
• Replace FastQuantities with DynamicQuantities.jl
• Tullio.jl/TensorCast.jl integration

Basic usage

• Add using UnitfulTensors at the top of your script
• Add units to your inputs and constants: ħ = 1.054571817e-34 * u"J*s"
• Create UnitfulTensors instead of arrays: UnitfulTensor([1u"s" 2u"m"; 3u"m" 4u"m^2/s"])
• Now linear algebra functions should just work

See Guide for more details.

Showcase

using UnitfulTensors, LinearAlgebra

# h-parameters
H = UnitfulTensor([3.7u"kΩ" 1.3e-4
                   120      8.7u"μS"])
Zₛ = 50u"Ω" # Source impedance
Yₗ = 20u"mS" # Load admittance
Hₛ = deepcopy(H)
Hₛ[1, 1] += Zₛ
Zout = inv(Hₛ)[2, 2] # Output impedance
Hₛₗ = H + Diagonal(UnitfulTensor([Zₛ, Yₗ]))
Gₜ = 4 * Zₛ * Yₗ * (Hₛₗ[2, 1] / det(Hₛₗ))^2 # Transducer power gain
Gₜ₂ = 4 * Zₛ * Yₗ * (inv(Hₛₗ)[2, 1])^2 # Equivalent expression
println((Zout, value(Zout), Zout/1u"kΩ"))
println((Gₜ, Gₜ₂))

# output

(220264.31718061675 kg m^2 A^-2 s^-3, 220264.31718061675, 220.26431718061676)
(10.2353526224919, 10.2353526224919)

Performance

Dense 3×3 matrix

using UnitfulTensors, BenchmarkTools
        
N = 3

A = randn(N, N)
b = randn(N)

uA = UnitfulTensor(A, 𝐋 * nodims(N, N))
ub = UnitfulTensor(b, 𝐌 * nodims(N))

println("Matrix * vector:")
@btime $A * $b
@btime $uA * $ub
println("Matrix * matrix:")
@btime $A * $A
@btime $uA * $uA
println("Solving a linear system:")
@btime $A \ $b
@btime $uA \ $ub

# output

Matrix * vector:
  119.563 ns (1 allocation: 80 bytes)
  160.231 ns (1 allocation: 80 bytes)
Matrix * matrix:
  82.297 ns (1 allocation: 128 bytes)
  115.884 ns (1 allocation: 128 bytes)
Solving a linear system:
  671.548 ns (3 allocations: 288 bytes)
  861.569 ns (5 allocations: 576 bytes)

Laplacian on a 10×10 grid, dense and sparse

using UnitfulTensors, LinearAlgebra, SparseArrays, BenchmarkTools

Nx = 10
N = Nx^2

Δ1d = Tridiagonal(fill(1., Nx - 1), fill(-2., Nx), fill(1., Nx - 1))
𝟙 = I(Nx)

Δ = kron(Δ1d, 𝟙) + kron(𝟙, Δ1d)
Δs = sparse(Δ)
ρ = randn(N)

uΔ = UnitfulTensor(Δ, 𝐋^-2 * nodims(N, N))
uΔs = UnitfulTensor(Δs, 𝐋^-2 * nodims(N, N))
uρ = UnitfulTensor(ρ * u"C/m^3")

println("Matrix * vector:")
@btime $Δ * $ρ
@btime $uΔ * $uρ
@btime $Δs * $ρ
@btime $uΔs * $uρ
println("Matrix * matrix:")
@btime $Δ * $Δ
@btime $uΔ * $uΔ
@btime $Δs * $Δs
@btime $uΔs * $uΔs
println("Solving a linear system:")
@btime $Δ \ $ρ
@btime $uΔ \ $uρ
@btime $Δs \ $ρ
@btime $uΔs \ $uρ

# output

Matrix * vector:
  1.618 μs (1 allocation: 896 bytes)
  2.038 μs (1 allocation: 896 bytes)
  755.763 ns (1 allocation: 896 bytes)
  1.128 μs (1 allocation: 896 bytes)
Matrix * matrix:
  56.211 μs (2 allocations: 78.17 KiB)
  57.104 μs (2 allocations: 78.17 KiB)
  16.128 μs (6 allocations: 35.50 KiB)
  16.492 μs (6 allocations: 35.50 KiB)
Solving a linear system:
  150.873 μs (4 allocations: 79.92 KiB)
  153.623 μs (6 allocations: 85.67 KiB)
  455.481 μs (125 allocations: 253.85 KiB)
  455.188 μs (127 allocations: 259.60 KiB)

Naively one could expect a 8x slowdown because 8 numbers are required to describe a single unitful scalar (1 for the numerical value, 7 for dimensions). This is not the case because

• The dimensions of a UnitfulTensor are of tensor product form, so only about 7*(M + N) numbers are required to describe the dimensions of an M×N matrix.
• Operations on dimensions are usually simpler then on numerical values. Multiplication of two N×N dense matrices requires O(N^3) operations on values and O(N) operations on dimensions.
• Most linear algebra functions reuse the dimensions of their arguments. For example, matrix multiplication copies row dimensions from the first argument and column dimensions from the last argument. In such cases, only the references to the corresponding AxisDimensions are copied from one UnitfulTensor to another, not the AxisDimensions themselves.
• Other operations invert AxisDimensions (inv, \) or introduce a dimensionally homogeneous axis (eigvecs), so they can be sped up by implementing NoAxisDimensions and InverseAxesDimensions (not implemented yet).
• Exponents of SI base units are usually rational numbers with small numerators and denominators, mostly integers and occasional half-integers. UnitfulTensors.jl stores them as Float32 instead of Float64, and even smaller options can be considered.
• Physical dimensions are usually far less diverse than numerical values, so memoization may be beneficial (not implemented yet).

The most popular Julia package for working with unitful quantities is Unitful.jl. It stores units in type parameters in an attempt to move unit computation to the compilation stage. This doesn't work well in the case of dimensionally heterogeneous arrays.

UnitfulTensors.jl has a FastQuantities submodule, which defines a type SIDimensions, which stores the physical dimensions of a UnitfulScalar as 7 Float32 numbers. The speedup is noticeable:

using UnitfulTensors, BenchmarkTools
using Unitful: @u_str as @uf_str

N = 1000

u = [1. * u"s"^i for i in 1:N]
v = [1. * uf"s"^i for i in 1:N] 

println("FastQuantities")
@btime $u .* $u
println("Unitful.jl")
@btime $v .* $v

# output

FastQuantities
  4.073 μs (2 allocations: 39.11 KiB)
Unitful.jl
  1.735 ms (1024 allocations: 32.44 KiB)

Brief API overview

UnitfulTensor represents an array of unitful quantities that can be used in linear or multilinear algebra. It stores an AbstractArray of numerical values and AxesDimensions that represent its physical dimensions.

For the numerical values, you can use basic Array, sparse matrices, StaticArrays, etc., but not OffsetArrays or anything else with non-standard indexing.

The physical dimensions must be of tensor product form, otherwise such an array couldn't be multiplied by any vector. For example, the dimensions of a UnitfulMatrix are the product of its row and column dimensions.

To ensure that this decomposition is unique, the row and column dimensions are normalized (divided by scalar factors so that the first row and column dimensions are NoDims). The combined normalization factor is stored as a separate field that can be accessed via dimscale, while the dimensions along each axis are stored as AxisDimensions objects and can be accessed via normdims.

Related work

• In Julia: Goretkin, Gebbie, more
• In C++: Withopf, Withopf, Withopf, Withopf, Withopf on Sea, Withopf, Withopf
• In the 90s: Hart