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

[ggebbie]

Thanks for the cool and promising package!

On Julia Discourse there seemed to be some interest in extending this package for "multidimensional analysis" (following George W. Hart) which I mostly interpret to be linear algebra utilities for the QuantityArray type. I don't know whether this fits with the vision of this package or whether it is best completed in a separate package. Any comments/advice on this PR are appreciated.

Here I make a start using the existing DynamicQuantitiesLinearAlgebraExt to add functionality for: inv, (\), svd, SVD, Diagonal, eigen, and det for QuantityArrays. All functions here use the "algebraic" interpretation of dimensional analysis.

Because QuantityArrays represent a dimensionally uniform matrix, the number of applicable linear algebra methods is fairly restricted. In the parlance of Multidimensional Analysis, an array would be any combination of values and dimensions/units in a tabular form. A matrix or multipliable matrix would be a restricted dimensional structure such that multiplication and other linear algebra functions are possible. In this light, the QuantityArray name for uniform matrices is somewhat confusing.

I think that other multipliable matrices (i.e., non-uniform) would best be represented by a different struct with the dimensional range and domain as two of the fields. This goes beyond an extension to this package and I plan to experiment on this idea in a separate package.

[github-actions[bot]]

Benchmark Results

	main	a7acc5a341e7ca...	main/a7acc5a341e7ca...
Quantity/creation/Quantity(x)	2.79 ± 0 ns	3.42 ± 0.011 ns	0.818
Quantity/creation/Quantity(x, length=y)	3.41 ± 0.01 ns	3.11 ± 0.01 ns	1.1
Quantity/with_numbers/*real	3.11 ± 0.01 ns	3.11 ± 0.01 ns	1
Quantity/with_numbers/^int	8.37 ± 2.5 ns	8.05 ± 2.2 ns	1.04
Quantity/with_numbers/^int * real	8.05 ± 2.2 ns	8.67 ± 2.5 ns	0.929
Quantity/with_quantity/+y	4.04 ± 0.01 ns	4.04 ± 0.001 ns	1
Quantity/with_quantity//y	3.11 ± 0.001 ns	3.42 ± 0.011 ns	0.909
Quantity/with_self/dimension	2.79 ± 0.001 ns	3.1 ± 0.01 ns	0.9
Quantity/with_self/inv	3.13 ± 0.01 ns	3.11 ± 0 ns	1.01
Quantity/with_self/ustrip	2.79 ± 0 ns	3.71 ± 0.92 ns	0.754
QuantityArray/broadcasting/multi_array_of_quantities	0.143 ± 0.00069 ms	0.144 ± 0.00094 ms	0.997
QuantityArray/broadcasting/multi_normal_array	0.0529 ± 0.0002 ms	0.0558 ± 0.0029 ms	0.947
QuantityArray/broadcasting/multi_quantity_array	0.157 ± 0.00076 ms	0.157 ± 0.00062 ms	1
QuantityArray/broadcasting/x^2_array_of_quantities	22.8 ± 2.1 μs	22.9 ± 2.3 μs	0.994
QuantityArray/broadcasting/x^2_normal_array	4.56 ± 1.2 μs	4.7 ± 1.2 μs	0.97
QuantityArray/broadcasting/x^2_quantity_array	6.95 ± 0.22 μs	6.96 ± 0.26 μs	0.999
QuantityArray/broadcasting/x^4_array_of_quantities	0.0844 ± 0.00043 ms	0.0814 ± 0.0006 ms	1.04
QuantityArray/broadcasting/x^4_normal_array	0.0498 ± 0.00016 ms	0.0497 ± 0.003 ms	1
QuantityArray/broadcasting/x^4_quantity_array	0.05 ± 0.00017 ms	0.053 ± 0.0029 ms	0.943
time_to_load	0.135 ± 0.0017 s	0.144 ± 0.0014 s	0.941

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR. Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[MilesCranmer]

Nice work!! I added some initial comments.

[ggebbie]

Failed Julia 1.6 runtest required some additional code for svd and Diagonal. I believe that the code should work for 1.9 and 1.6 (awaiting 1.6 GitHub Action to be sure). Involved adding some uglier code, however.

[MilesCranmer]

Edit: I guess we might need to add a special method for diagm as well, since it uses zero(::Type{T}) rather than zero(::T). So best to ustrip -> diagm -> QuantityArray(result, dimension(input))

[ggebbie]

Edit: I guess we might need to add a special method for diagm as well, since it uses zero(::Type{T}) rather than zero(::T). So best to ustrip -> diagm -> QuantityArray(result, dimension(input))

This is the same issue that led me to writing Diagonal. Right -- best to make a special diagm function.

[MilesCranmer]

I guess some of these issues could also be fixed if we could figure out a way to implement https://github.com/SymbolicML/DynamicQuantities.jl/issues/76, right? Trading type instability with compatibility (though people could get type stability if they simply used zero(::T) instead... but at least it would avoid any errors).

[ggebbie]

Getting closer and hope to not have missed anything. Certainly ok if you edit this PR. I tried to test on Julia 1.6 locally but Pkg complained about Compat... package not being in registry so I quit.

[ggebbie]

Compatibility/composability with other packages is a big reason why I'm interested in this package. Unitful doesn't work with Turing, for example. It seems like there are issues with SciML and Unitful as well. Making some design choices in this package so that it is more usable in combination with other packages would be a good idea in my opinion.

[MilesCranmer]

Hey @ggebbie,

I just made a few changes:

1. I incorporated the new abstract interface to your LinearAlgebra functions, so that any AbstractQuantity or AbstractGenericQuantity is permitted, rather than just Quantity. (Uses quantity_type to get the base quantity out of QuantityArray)
2. I removed the where T statements to simplify things. I realized that those errors will just get caught anyways by the method operating on the raw array, so there's no need to limit it there.
3. Switched the dimension errors to use DimensionError rather than ErrorException.
4. I added Base.:/.
5. I introduced promotion of quantity types in *, /, \. This means you can multiply an array of symbolic units with one of SI units, like :

julia> x = QuantityArray(rand(3, 3), u"km/s"); typeof(x)
QuantityArray(::Matrix{Float64}, ::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}})

julia> y = QuantityArray(rand(3, 3), us"km/s"); typeof(y)
QuantityArray(::Matrix{Float64}, ::Quantity{Float64, SymbolicDimensions{DynamicQuantities.FixedRational{Int32, 25200}}})

julia> x * y |> typeof
QuantityArray(::Matrix{Float64}, ::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}})

Let me know if this sounds good. Also, I have a question for you: is the output dimension of eigen supposed to be dimension(A)? Shouldn't the output of eigen just be a regular scalar for the eigenvalue, and a regular vector for the eigenvectors? Or maybe it is good the way it is?? Just wanted to double check.

If this is true then eigen should be changed to reflect this.

Also, could you add unittests for the remaining pieces? I aim to have 100% coverage of all the code in the repo. (As well as checking that errors are caught with @test_throws)

Thanks! Miles

[MilesCranmer]

@j-fu this is probably also of interest to you

[ggebbie]

Looks like some good additions above, but it will take me about a week to take a closer look due to holidays.

I might be misunderstanding your thought process with regards to eigen. I checked the code and it looks correct to me and is consistent with Hart's textbook. Consider the eigenvalue equation: Ax = lambda x. If A is uniform (i.e., all entries have the same dimension/unit), then the action of A increments (i.e., multiplies) the dimension of x by dimension(A). For the right hand side of the equation to have a consistent dimension with the left hand side, dimension(lambda) must be equal to dimension(A). x could be any uniform vector, but the convention is to make it dimensionless, as it is the simplest eigenvector that satisfies the equation. So it makes sense to report the eigenvalues with units but the eigenvectors without.

In more complicated (non-uniform) matrices, both lambda and x may require units. And for other non-uniform matrices, there may be no eigenstructure at all. The requirement is that a matrix is squarable (i.e., A*A is dimensionally conformable) for eigenvalues and eigenvectors to exist.

[MilesCranmer]

Sounds good! Thanks for clarifying

[MilesCranmer]

Ping on this @ggebbie; what's the status and are there any blockers I can help with?

[ggebbie]

Glad that this is still on the radar. I haven't found time to work on it. It looks like the action items are 1) improving test coverage and 2) implementing @test_throws. I can certainly handle 1) and will let you know if 2) is a blocker after I learn more about it.

[BambOoxX]

Just for reference. https://discourse.julialang.org/t/ann-unitfultensors-jl-efficient-arrays-with-units/114028/4

[MilesCranmer]

@ggebbie, do you have time to finish it off? It seems like you were almost done. Even if not all of the interface is ready, that would be much appreciated. Like inv and \ would be great!