Compatibility with backslash and sparse arrays

[briochemc]

So I am not sure what the issue is, or if I am not doing things properly, but it seems that it should be doable to make Unitful.jl work with linear algebra. I tend to do array multiplications and solve linear systems a lot via backslash. Here is an example of unit-less things I do:

# no-unit Array example that works
κ = rand(3, 3)
c = rand(3)
s = κ * c
c ≈ κ \ s

# no-unit sparse Array example that works
using SparseArrays, LinearAlgebra
κ = sprand(3, 3, 0.5) + I
c = rand(3)
s = κ * c
c ≈ κ \ s

Now doing that with unit-full scalars is not an issue:

using Unitful
# Scalar example that works
κ = rand()u"1/s"
c = rand()u"mol/m^3"
s = κ * c
c ≈ κ \ s

But it breaks on two important things for me:

1. Using \ on arrays raises an error:

   # Array example that works for multiplication only
   κ = rand(3, 3)u"1/s"
   c = rand(3 )u"mol/m^3"
   s = κ * c
   c ≈ κ \ s # Gives DimensionError

   My guess is that the factorization is a problem, as I do not know how would one split the units of κ into the units of its factors. (E.g., what are the units of L and U in κ = L * U?). Is it possible to give a single unit to a whole array and tell \ to do the unit conversion on the side? (I.e., do the κ \ s evaluation without units, and assign it the unit of s divided by the unit of κ.)

2. Using * on sparse arrays raises an error:

   # Sparse Array example that does not work
   κ = (sprand(3, 3, 0.5) + I)u"1/s"
   c = rand(3)u"mol/m^3"
   s = κ * c # Gives DimensionError

   My guess here is it comes from the zeros of the sparse matrix, which are not given a unit? Maybe doing the unit conversion on the side is still doable?

EDIT: I realize this is a partial duplicate of issue #46. The only new feature request here is thus the sparse matrix one.

[briochemc]

So I think I started a solution, based on @ChrisRackauckas suggestion. But it is ugly so I'm sure someone could help clean this up, and maybe finish the work with factorize? Hopefully this is useful 😄

using Unitful, LinearAlgebra, SparseArrays

# Overload \ to work with unitful matrices
function Base.:\(A::Array{Quantity{TA,DA,UA}, 2}, x::Vector{Tx}) where {TA,DA,UA,Tx}
     return (ustrip.(A) \ ustrip.(x)) * (unit(x[1]) / unit(A[1]))
end

# Overload * to work with either or both A and x unitful (and A is sparse)
function Base.:*(A::SparseMatrixCSC{Quantity{TA,DA,UA}, Int64}, x::Vector{Quantity{Tx,Dx,Ux}}) where {TA,DA,UA,Tx,Dx,Ux}
     return (ustrip.(A) * ustrip.(x)) * (unit(A.nzval[1]) * unit(x[1]))
end
function Base.:*(A::SparseMatrixCSC{TA, Int64}, x::Vector{Quantity{Tx,Dx,Ux}}) where {TA,Tx,Dx,Ux}
     return (ustrip.(A) * ustrip.(x)) * (unit(A.nzval[1]) * unit(x[1]))
end
function Base.:*(A::SparseMatrixCSC{Quantity{TA,DA,UA}, Int64}, x::Vector{Tx}) where {TA,DA,UA,Tx}
     return (ustrip.(A) * ustrip.(x)) * (unit(A.nzval[1]) * unit(x[1]))
end

# Overload \ to work with unitful A or x (or both) and A is sparse
function Base.:\(A::SparseMatrixCSC{Quantity{TA,DA,UA}, Int64}, x::Vector{Quantity{Tx,Dx,Ux}}) where {TA,DA,UA,Tx,Dx,Ux}
     return (ustrip.(A) \ ustrip.(x)) * (unit(x[1]) / unit(A.nzval[1]))
end
function Base.:\(A::SparseMatrixCSC{Quantity{TA,DA,UA}, Int64}, x::Vector{Tx}) where {TA,DA,UA,Tx}
     return (ustrip.(A) \ ustrip.(x)) * (unit(x[1]) / unit(A.nzval[1]))
end
function Base.:\(A::SparseMatrixCSC{TA, Int64}, x::Vector{Quantity{Tx,Dx,Ux}}) where {TA,Tx,Dx,Ux}
     return (ustrip.(A) \ ustrip.(x)) * (unit(x[1]) / unit(A.nzval[1]))
end

# ---------- Non sparse tests ----------
A = rand(3, 3) + I
Au = A * 1u"1/s"
x = rand(3)
xu = x * 1u"mol/m^3"
# Test *
A * x
Au * xu
Au * x
A * xu
# Test \
A \ x
Au \ x
A \ xu
Au \ xu
# ---------- Sparse tests ----------
A = sprand(3, 3, 0.5) + I
Au = A * 1u"1/s"
x = rand(3)
xu = x * 1u"mol/m^3"
# Test *
A * x
Au * x
A * xu
Au * xu
# Test \
A \ x
Au \ x
A \ xu
Au \ xu

One thing missing (that I personally need) would be to overload factorize to only assign the unit to, e.g., the first factor, so that then \ would work with factorize(Au) \ xu and other combinations, I guess...

[ggebbie]

This comment is more than 4.5 years old, but I see some activity in 2022, so I've given a shot at addressing the issue. I've written a package at https://github.com/ggebbie/UnitfulLinearAlgebra.jl that can pass the tests with a couple of added lines. Thanks to @briochemc for getting this kickstarted with his code above.

• backslash should generally work
• LU factorization gives an object with the right units
• sparse matrices are partially handled, but not in the way described above
• unitful identity matrices are constructable, but require more information

Any suggestions/help on the package are appreciated, including a proper name.

Here are the tests from above with comments:

ENV["UNITFUL_FANCY_EXPONENTS"] = true
using Revise
using UnitfulLinearAlgebra
using Unitful
using LinearAlgebra
using Test
            A = rand(3, 3) + I
            Au = A * 1u"1/s"

            # A with multipliable matrix
            Amm = BestMultipliableMatrix(Au)

            x = rand(3)
            xu = x * 1u"mol/m^3"
            # Test *
            A * x
            Au * xu
            Au * x
            A * xu
            # Test \
            A \ x # works with a UniformMatrix or LeftUnitformMatrix
            #Au \ x # won't work
            Amm \ x # gets units right
            #A \ xu # won't work
            #Au \ xu # no existing method
            Amm \ xu

            # ---------- Sparse tests ----------
            A = sprand(3, 3, 0.5) + I
            Au = A * 1u"1/s"
            Ammfull = BestMultipliableMatrix(Matrix(Au))# not working with SparseArray now
            Amm = BestMultipliableMatrix(A,fill(u"mol/m^3",3),fill(u"s*mol/m^3",3))  # use constructor, internally stores a sparse matrix
            x = rand(3)
            xu = x * 1u"mol/m^3"

            # Test *
            A * x
            Au * x
            A * xu
            Au * xu
            Amm* xu
            # Test \

            # Problem: units not right for x to be conformable with Au.
            # change x to y
            y = rand(3);
            yu = y.*unitrange(Amm)
            A \ y 
            #Au \ x # stack overflow, doesn't work at lu, no method
            Amm \ y # is UniformMatrix, so it works
            #A \ yu # doesn't work, no method
            #Au \ yu, doens't work, no lu method
            Amm \ yu # works, same units as x

[goretkin]

Just want to point out some related work: https://github.com/goretkin/UnitfulLinearAlgebra.jl

It's a (immature) approach where unit information is encoded in the type itself. This is probably preferable when the arrays are small and backed by arrays in StaticArrays.jl

On Wed, Dec 7, 2022, 11:51 AM G Jake Gebbie @.***> wrote:

This comment is more than 4.5 years old, but I see some activity in 2022, so I've given a shot at addressing the issue. I've written a package at https://github.com/ggebbie/UnitfulLinearAlgebra.jl that can pass the tests with a couple of added lines. Thanks to @briochemc https://github.com/briochemc for getting this kickstarted with his code above.

• backslash should generally work
• LU factorization gives an object with the right units
• sparse matrices are partially handled, but not in the way described above
• unitful identity matrices are constructable, but require more information

Any suggestions/help on the package are appreciated, including a proper name.

Here are the tests from above with comments:

ENV["UNITFUL_FANCY_EXPONENTS"] = true using Revise using UnitfulLinearAlgebra using Unitful using LinearAlgebra using Test A = rand(3, 3) + I Au = A * 1u"1/s"

        # A with multipliable matrix
        Amm = BestMultipliableMatrix(Au)

        x = rand(3)
        xu = x * 1u"mol/m^3"
        # Test *
        A * x
        Au * xu
        Au * x
        A * xu
        # Test \
        A \ x # works with a UniformMatrix or LeftUnitformMatrix
        #Au \ x # won't work
        Amm \ x # gets units right
        #A \ xu # won't work
        #Au \ xu # no existing method
        Amm \ xu

        # ---------- Sparse tests ----------
        A = sprand(3, 3, 0.5) + I
        Au = A * 1u"1/s"
        Ammfull = BestMultipliableMatrix(Matrix(Au))# not working with SparseArray now
        Amm = BestMultipliableMatrix(A,fill(u"mol/m^3",3),fill(u"s*mol/m^3",3))  # use constructor, internally stores a sparse matrix
        x = rand(3)
        xu = x * 1u"mol/m^3"

        # Test *
        A * x
        Au * x
        A * xu
        Au * xu
        Amm* xu
        # Test \

        # Problem: units not right for x to be conformable with Au.
        # change x to y
        y = rand(3);
        yu = y.*unitrange(Amm)
        A \ y
        #Au \ x # stack overflow, doesn't work at lu, no method
        Amm \ y # is UniformMatrix, so it works
        #A \ yu # doesn't work, no method
        #Au \ yu, doens't work, no lu method
        Amm \ yu # works, same units as x

— Reply to this email directly, view it on GitHub https://github.com/PainterQubits/Unitful.jl/issues/150#issuecomment-1341266936, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEN3M7732XOCBGZHRYNTQTWMC6A3ANCNFSM4FSILVBA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[ggebbie]

thanks a lot @goretkin for posting your work. Thanks also for the references in your README which were helpful for my learning. I have put the unit information into attributes of the MultipliableMatrix struct. Only the dimensional domain and range are needed, which should help with performance. At this point, however, it looks like I am incurring some performance penalty when doing matrix multiplication. Factorizations don't seem to be too much slower than the unitless version, on the other hand.

[goretkin]

I think I prefer your approach of having fields corresponding to unit information. This does not necessarily have to incur a memory or performance penalty in the small array case by using singleton values, e.g. Using Val{8}() instead of 8.

It will require experimentation, but I think this is the most flexible design.

On Wed, Dec 7, 2022, 2:12 PM G Jake Gebbie @.***> wrote:

thanks a lot @goretkin https://github.com/goretkin for posting your work. Thanks also for the references in your README which were helpful for my learning. I have put the unit information into attributes of the MultipliableMatrix struct. Only the dimensional domain and range are needed, which should help with performance. At this point, however, it looks like I am incurring some performance penalty when doing matrix multiplication. Factorizations don't seem to be too much slower than the unitless version, on the other hand.

— Reply to this email directly, view it on GitHub https://github.com/PainterQubits/Unitful.jl/issues/150#issuecomment-1341464750, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAEN3M464JDQERGPIRKZH6LWMDOSRANCNFSM4FSILVBA . You are receiving this because you were mentioned.Message ID: @.***>