Linear Algebra extension: svd, inv, *, /, \ for QuantityArray

ext/DynamicQuantitiesLinearAlgebraExt.jl

@@ -1,8 +1,69 @@
 module DynamicQuantitiesLinearAlgebraExt
 
-import LinearAlgebra: norm
-import DynamicQuantities: AbstractQuantity, ustrip, dimension, new_quantity
+import LinearAlgebra: norm, inv, (\), svd, Algorithm, default_svd_alg, SVD, Diagonal, Adjoint, Transpose, AbstractRotation, AbstractMatrix, eigen, eigsortby, Eigen, det
+import DynamicQuantities: AbstractQuantity, ustrip, dimension, new_quantity, AbstractDimensions, QuantityArray, Quantity
+
+const QuantityArrayVecOrMat{T} = Union{QuantityArray{T,2},QuantityArray{T,1}} where T
 
 norm(q::AbstractQuantity, p::Real=2) = new_quantity(typeof(q), norm(ustrip(q), p), dimension(q))
 
+\(q::QuantityArrayVecOrMat,r::QuantityArrayVecOrMat) = QuantityArray(ustrip(q)\ustrip(r),dimension(r)/dimension(q))
+\(q::QuantityArrayVecOrMat,r::Union{AbstractVector,AbstractMatrix}) = QuantityArray(ustrip(q)\r,inv(dimension(q)))
+# not implemented, AbstractMatrix \ QuantityArray
+
+inv(Q::QuantityArray{T,2}) where T = QuantityArray(inv(ustrip(Q)),inv(dimension(Q)))
+
+"""
+    svd(A::QuantityArray; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD
+
+    Singular value decomposition (SVD) of `QuantityArray`.
+    Exists for uniform matrices which includes all `QuantityArray`s (pp. 124, Hart, 1995).
+
+    Returns SVD factorization of parametric type: `SVD{T, Quantity{T, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Matrix{T}, QuantityArray{T, 1, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}, Quantity{T, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Vector{T}}}`.
+
+    Factorization `F` can be deconstructed: `U,σ,V = F`. 
+
+    Functions working: , `size`, `adjoint`.
+    Functions partially working: `inv`
+    Functions not working: `svdvals`, `ldiv!`.
+"""
+function svd(A::QuantityArray;full=false,alg::Algorithm = default_svd_alg(ustrip(A))) 
+    F = svd(ustrip(A), full=full, alg=alg)
+    #return SVD(F.U,QuantityArray(F.S,dimension(A)),F.Vt) # julia 1.9 pass but 1.6 fail because 2nd argument is not a Vector
+    S = [Quantity(F.S[i],dimension(A)) for i in eachindex(F.S)]
+    return SVD(F.U,S,F.Vt) # julia 1.6 pass but long-winded
+end
+
+Diagonal(q::QuantityArray)  = QuantityArray(Diagonal(ustrip(q)),dimension(q))
+function Diagonal(q::Vector{Quantity{T,D}}) where {T,R,D<:AbstractDimensions{R}}
+    if allequal(dimension.(q))
+        return QuantityArray(Diagonal(ustrip.(q)),dimension(q[begin]))
+    else
+        error("Diagonal matrix would not be a `QuantityArray`. Vector must contain `Quantity`s with same dimension.")
+    end
+end
+
+"""
+    function eigen(A::T;permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T = QuantityArray
+
+    Thin wrapper for `eigen` with same keyword arguments as `LinearAlgebra.eigen`.
+    Only squarable matrices have eigenstructure (pp. 96, Hart, 1995).
+    Eigenvalues have the same dimensions as 𝐀[1,1].
+    Eigenvectors are parallel to the domain and range.
+    There are multiple ways to distribute the units among the eigenvectors, however.
+    If 𝐀 is endomorphic (i.e., the dimensional domain and range are the same), then the dimensional domain should
+    be taken as the units of the eigenvectors (pp. 205, Hart, 1995).  
+    In the general case, physical intuition and the equation 𝐀𝐱 = λ𝐱
+    dictate that the units of the eigenvectors are equal to the dimensional domain of 𝐀 (pp. 206, Hart, 1995).
+
+    The following functions are available for `Eigen(::QuantityArray)` objects: eigvals, [`det`](@ref).  
+    Functions not working: [`inv`](@ref) and [`isposdef`](@ref).
+"""
+function eigen(A::QuantityArray;permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) 
+    F = eigen(ustrip(A), permute=permute, scale=scale, sortby=sortby)
+    return Eigen(QuantityArray(F.values,dimension(A)), F.vectors)
+end
+
+det(A::QuantityArray) = Quantity(det(ustrip(A)),dimension(A)^(size(A,1)))
+
 end

src/arrays.jl

@@ -212,3 +212,10 @@ ucurrent(q::QuantityArray) = ucurrent(dimension(q))
 utemperature(q::QuantityArray) = utemperature(dimension(q))
 uluminosity(q::QuantityArray) = uluminosity(dimension(q))
 uamount(q::QuantityArray) = uamount(dimension(q))
+
+# handle ambiguities in multiplication
+Base.:*(q::QuantityArray,r::QuantityArray) = QuantityArray(ustrip(q)*ustrip(r),dimension(q)*dimension(r))
+Base.:*(q::QuantityArray,r::AbstractMatrix) = QuantityArray(ustrip(q)*r,dimension(q))
+Base.:*(q::QuantityArray,r::AbstractVector)  = QuantityArray(ustrip(q)*r,dimension(q))
+Base.:*(r::AbstractMatrix,q::QuantityArray) = QuantityArray(r*ustrip(q),dimension(q))
+Base.:*(r::AbstractVector,q::QuantityArray) = QuantityArray(r*ustrip(q),dimension(q))

test/runtests.jl

@@ -22,6 +22,9 @@ else
     @safetestset "Unit tests" begin
         include("unittests.jl")
     end
+    @safetestset "LinearAlgebra.jl integration tests" begin
+        include("test_linearalgebra.jl")
+    end
     @safetestset "Aqua tests" begin
         include("test_aqua.jl")
     end

test/test_linearalgebra.jl

@@ -0,0 +1,78 @@
+#import Pkg; Pkg.activate(".")
+#using Revise
+using DynamicQuantities, LinearAlgebra
+
+# checks for unit consistency with subtraction
+within(A,B,tol) =  maximum(abs.(ustrip(A - B))) < tol
+
+v = randn(2,2)
+d = 1u"m"
+e = u"m"
+A = QuantityArray(v, d) # Create a `QuantityArray` with value `value` and dimensions `dimensions`.
+
+# vector multiplication
+q = [1,2]
+r = A\q
+@test isequal(dimension(A*r),dimension(QuantityArray(q)))
+@test within(ustrip(A*r),q,1e-10)
+
+# ldiv 2D QuantityArray and 1D QuantityArray
+g = QuantityArray(randn(2), d) 
+h = A\g
+@test within(A*h,g,1e-10)
+
+# test inv
+B = inv(A)
+@test isapprox(B*A,I(2),rtol = 1e-10)
+@test isapprox(A*B,I(2),rtol = 1e-10)
+
+vc = randn(2,2)
+dc = 1u"kg"
+C = QuantityArray(vc, dc) # Create a `QuantityArray` with v
+
+# check ldiv with 2D x 2D `QuantityArray`s
+D = A*C
+@test within(A\D,C,1e-12)
+
+F = svd(A)
+U,σ,V = F
+
+inv(F) # doesn't return a QuantityArray, could add a wrapper here
+
+@test isequal(dimension(σ),dimension(A))
+
+Σ = Diagonal(σ)
+@test within(U*Σ*V',A,1e-10)
+
+# SVD in modal form
+# truncated SVD
+K = 2
+U[:,1:K]*Diagonal(σ[1:K])*V[:,1:K]'
+
+K = 1
+U[:,1:K]*Diagonal(σ[1:K])*V[:,1:K]'
+
+# retain QuantityArray, best way to store this info
+@test Σ*Σ isa QuantityArray
+@test U*Σ isa QuantityArray
+@test U'*Σ isa QuantityArray
+@test transpose(U)*Σ isa QuantityArray
+@test Σ*V isa QuantityArray
+@test Σ*V' isa QuantityArray
+@test Σ*transpose(V) isa QuantityArray
+@test U*Σ*V' isa QuantityArray
+
+F = eigen(A)
+vals,vecs = F
+@test isequal(dimension(vals),dimension(A))
+# inv(F) # doesn't work
+@test within(det(F),det(A),1e-10)
+
+# eigenstructure satisfied?
+@test within(A*vecs,vecs*Diagonal(vals),1e-10)
+
+det(A) # does it run?
+
+# svdvals(A) # error, LinearAlgebra.eigen.jl doesn't handle it
+
+# isapprox(ustrip(transpose(A)),transpose(ustrip(A))) # no show method for Transpose{QuantityArray}