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 Base: (*) #, transpose
+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
 
 norm(q::AbstractQuantity, p::Real=2) = new_quantity(typeof(q), norm(ustrip(q), p), dimension(q))
 
+# handle ambiguities
+*(q::QuantityArray,r::QuantityArray) = QuantityArray(ustrip(q)*ustrip(r),dimension(q)*dimension(r))
+*(q::QuantityArray,r::AbstractMatrix) = QuantityArray(ustrip(q)*r,dimension(q))
+*(q::QuantityArray,r::AbstractVector)  = QuantityArray(ustrip(q)*r,dimension(q))
+*(r::AbstractMatrix,q::QuantityArray) = QuantityArray(r*ustrip(q),dimension(q))
+*(r::AbstractVector,q::QuantityArray) = QuantityArray(r*ustrip(q),dimension(q))
+
+# a choice to make this function type stable (to be reviewed/revisited)
+#transpose(q::QuantityArray) = QuantityArray(transpose(ustrip(q)),dimension(q))
+
+\(q::QuantityArray,r::QuantityArray) = QuantityArray(ustrip(q)\ustrip(r),dimension(r)/dimension(q))
+\(q::QuantityArray,r::Union{AbstractVector,AbstractMatrix}) = QuantityArray(ustrip(q)\r,dimension(q)^-1)
+# not implemented, AbstractMatrix \ QuantityArray
+
+inv(Q::QuantityArray) = QuantityArray(inv(ustrip(Q)),dimension(Q)^-1)
+
+"""
+    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)
+end
+
+Diagonal(q::QuantityArray)  = QuantityArray(Diagonal(ustrip(q)),dimension(q))
+
+"""
+    function eigen(A::T;permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T <: AbstractMultipliableMatrix
+
+    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 A[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

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,74 @@
+#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)
+
+# 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
+
+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, eigen.jl doesn't handle it
+
+# isapprox(ustrip(transpose(A)),transpose(ustrip(A))) # no show method for Transpose{QuantityArray}