JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.85k stars 5.49k forks source link

Feature request: Decompositions for Unitful matrices #56659

Open NAThompson opened 1 day ago

NAThompson commented 1 day ago

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 commented 1 day ago

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