Is there a \ operator for Unitful?

I was running into a problem using the \ (linear system solver) operator in Julia. I guess I wrongly assumed the \ operator was overloaded for matrices/vectors with entries of type <: Number (witch unitful types are). Is there any way to easily overload the \ operator to accept unitful types, or perhaps a different operator/function that works for all types that are subtypes of Number (since I plan on using Sym and Float/Integer types in my code as well)

Environment: Pluto on MacOS
Packages: Unitful, Sympy, Plots, PlutoUI

	a11 = 1u"m"
	a12 = 2u"ft"
	a21 = 5u"mm"
	a22 = 0u"inch"
	A = [a11 a12; a21 a22]
	b = [3u"inch"; 1u"mm"]
	x  = A \ b # This line throws error

MethodError: no method matching (Unitful.Quantity{Rational{Int64}})(::Rational{Int64})

1 Like

I don’t think there is an easy way to define a unique, unambiguous solution to that example problem. Which should the units of the elements of x? Moreover, in this case it’s only a matter of scale between units of the same kind, but for arbitrary units in A and b, it might be impossible to solve.

1 Like

Maybe something along these lines:

using LinearAlgebra, Unitful
import LinearAlgebra: \

(\)(A::Matrix{Quantity{T,U,V}}, B::AbstractVecOrMat) where {T,U,V} =
  Quantity.(ustrip(A)\ustrip(B), unit(eltype(B))/unit(eltype(A)))

With this definition:

julia> A\b
2-element Vector{Rational{Int64}}:

julia> c = ustrip(b)
2-element reinterpret(Rational{Int64}...):

julia> A\c
2-element Vector{Quantity{Rational{Int64}...}}:
       1//5 m^-1
 -619//3048 m^-1

But this was a quick sketch, might be ins-and-outs to fix yet.

Oh yeah, it makes sense that the ‘b’ vector would be in units of Length*mysteryunit since Ax = b, and x is unknown. Maybe x can be solved given the units of b?

Using the one-argument ustrip method is generally a bad idea: you don’t know what you’re stripping, it works well only for the trivial case where you have all homogeneous quantities.


Thanks for the advice! I changed up the method a little bit to ensure that all unit of elements in the matrices A and B are the same before using the system solver.

function Base.:\(A::Matrix{Quantity{T,U,V}}, B::AbstractVecOrMat) where {T,U,V}
		if typeof(unit(eltype(A))) ≠ Unitful.FreeUnits{(), NoDims, nothing}
			A = map(x -> uconvert(unit(eltype(A)), x), A) 
			B = map(x -> uconvert(unit(eltype(B)), x), B)
			return Quantity.(ustrip(A)\ustrip(B), unit(eltype(B))/unit(eltype(A)))

There should be another definition of \ covering the case when A is a non-Unitful matrix and B is Unitful. Should be pretty similar to the one already defined.