Feature request: Decompositions for Unitful matrices

[NAThompson]

To reproduce:

julia> using Unitful
julia> using LinearAlgebra
julia> matrix = rand(4, 3) .* u"kg"
4×3 Matrix{Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}}:
 0.903809 kg  0.857622 kg  0.815957 kg
 0.959174 kg  0.235948 kg  0.379886 kg
 0.123471 kg   0.72895 kg  0.551869 kg
 0.573551 kg  0.819126 kg  0.911172 kg
julia> svd(matrix)
ERROR: MethodError: no method matching svd!(::Matrix{Quantity{Float64}}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
The function `svd!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  svd!(::Bidiagonal{var"#s4971", V} where {var"#s4971"<:Union{Float32, Float64}, V<:AbstractVector{var"#s4971"}}; full) got unsupported keyword argument "alg"
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/bidiag.jl:254
  svd!(::StridedMatrix{T}, ::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF64, ComplexF32} got unsupported keyword arguments "full", "alg"
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/svd.jl:395
  svd!(::StridedVector{T}; full, alg) where T<:Union{Float32, Float64, ComplexF64, ComplexF32}
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/svd.jl:109
  ...

Stacktrace:
 [1] svd(A::Matrix{Quantity{Float64, 𝐌, Unitful.FreeUnits{…}}}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/svd.jl:179
 [2] svd(A::Matrix{Quantity{Float64, 𝐌, Unitful.FreeUnits{(kg,), 𝐌, nothing}}})
   @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.11/LinearAlgebra/src/svd.jl:178
 [3] top-level scope
   @ REPL[11]:1
Some type information was truncated. Use `show(err)` to see complete types.

It appears that there are a few methods in the LinearAlgebra package which are fine with Unitful matrices (like Diagonal, transpose, Adjoint), but LU, Cholesky, and svd do not support it.

[giordano]

For background, there's a lengthy discussion in https://github.com/PainterQubits/Unitful.jl/issues/46

[timholy]

https://github.com/anonymous-shrew/UnitfulTensors.jl. I haven't used it, but it looks like it has the right design. You can't do matrix algebra operations unless a matrix's units are of the form u ⊗ v (i.e., at most a rank-1 matrix), as otherwise even matrix multiplication cannot be defined. So it really doesn't make much sense to try to implement generic algorithms for matrices of unrestricted unit combinations that will be slow and are not guaranteed to succeed. The better option is to coerce the inputs to "validated" types, which is what that package appears to do.

[aplavin]

Concretely typed Unitful matrices/vectors have the same units for each element. This is a common case that would be nice to support, and it doesn't require any of the more involved packages. AFAIU, UnitfulTensors.jl is only needed when one needs different units for rows/columns.

[andreasnoack]

See also https://github.com/JuliaLang/LinearAlgebra.jl/issues/128 and the linked issues

[NAThompson]

@timholy : Just to give additional context: This ask arose from an optimization routine where the matrix construction and decomposition is an implementation detail; the objective function F(x) can have dimensions in both arguments and output yielding a dimensioned matrix. Obviously I could request the user use UnitfulTensors.jl, but your package is the most popular, and it does seem reasonable to see how far the generic code can be pushed before making a user request.

Personally, I also find writing the optimization routines easier if I can test it with dimensioned quantities-speed be damned if need be.