handling of vector quantities and linear algebra libraries

[rconde01]

I'm curious to your thoughts on how this applies beyond scalar quantities. How might this be combined with a linear algebra library like eigen?

[mpusz]

Hi, I do not have much experience with Eigen but I can imagine 2 approaches here:

1. quantity behaves like a numeric type so it can be easily used in other libraries as a vector of quantities.
2. quantity takes Rep as a template parameter and maybe it could be some vector representation there.

I asked @hatcat if he could take a look at it. Maybe we should cooperate while working on units and LA libraries...

[hatcat]

Sorry about the delay. I have an insane three weeks coming up (C++ Russia in St Petersburg, WG21 in Belfast, Meeting C++ in Berlin).

Looking through the quantity class I see a lot of LA-like stuff going on (although I am confused by operator/). I think both approaches are appropriate: would you like to schedule a telecon to refine the discussion about which would be more appropriate under which conditions?

[mpusz]

As it is only 2 weeks to Belfast now, maybe we should meet there and discuss possible approaches here?

[mpusz]

I met with linear algebra folks in Belfast and the outcome is that the current LA proposal will work great with units library :-)

The only thing missing in LA is matrix_cast<To>(From) and an extension point for casting which will be called by it and will allow calling units::quantity_cast<To>(from) for each element of the matrix.

We agreed that we will provide examples of how both libraries can cooperate with each other.

[rconde01]

Great to hear. I'm very interested in seeing some examples. Some thoughts I have at the moment:

• IIRC the LA proposal is about the base operation implementation. It would be good to know it also works for Expression-Template layers (e.g. Eigen, Blaze) which presumably could layer over that to avoid temporaries. I have no particular reason to think it wouldn't work in this context, but I appreciate the possibility that there might be a blocker there.
• did your discussions include mixed units or only uniform units?
• maybe it's straightforward, but i wonder about unit erasing operations like diagonailization
• I suppose something like a rotation would be a unitless transformation?

Thanks for your work!

[dwith-ts]

@rconde01 : A rotation is a linear transformation between vector spaces and as such is very similar to a jacobian matrix (actually the jacobian of an isometry transformation is the rotation matrix).
=> rotations have unit-less entries

[dwith-ts]

Regarding the idea to put physical units inside a LA library:
This is a godd idea, but I presume that this can break pretty easily because all matrix operations would need to be supported.

This might be easy (if the linalg library is designed very specifically for this) for something like the dot product of 2 vectors where the return type is just decltype(vecA(0) * vecB(0)).

If you have matrices, there are very few matrices where you can get along by just specifying one element unit. A covariance matrix with unit [m^2] for example is something completely different than a Jacobian matrix with unit [m^2] and the allowed operations are also different for a covariance vs. a jacobian.

@rconde01 's question with a rotation is again another level where you would ideally also need support for specifying coordinate systems, not only for si-units.

Furthermore, the distinction between points and displacement vector is important: for example, a position vector can only be transformed with an isometry (rotation + translation) whereas a displacement vector can only be transformed with a rotation (no translation occurs).

To summarize, IMHO it is very valuable to represent units in linalg vectors and matrices, but this needs a unified treatment to also support non-uniform units in vectors and matrices and is best done by an additional layer on top of a LA library that has a user-interface which uses units but internally uses plain scalar vectors and matrices from any linalg library.

This also makes the LA core code easier on the compilers, a library like e.g. Eigen already contains lots of templates.

[kwikius]

Here FWIW is an example of a matrix multiply with physical quantities in my library

https://github.com/kwikius/quan-trunk/blob/master/quan_matters/examples/fusion/kalman1.cpp
(The matrix is essentially a tuple and, being a tuple, the compiler must precompute the type of the resulting matrix :
https://github.com/kwikius/quan-trunk/blob/master/quan/fun/make_matrix_mux_result.hpp#L140

I have done multiplication on 4,4 x 4,4 matrices and have tried larger matrices but compile time got very long and large amounts of memory where used by the compiler as the matrices got bigger

I also experimented with a "static value" which was designed as a matrix element to optimise out operations on constants by turning them into different types
https://github.com/kwikius/quan-trunk/tree/master/quan/fusion/static_value

An example optimisation is static_value<T1,0> * runtime_x -> static_value<decltype(T1{}*runtime_x),0> which is a nop at runtime and will be computed in the matrix_mux_result metafunction above.

Also compiletime /constexpr operations on such entities e.g inner_product
https://github.com/kwikius/quan-trunk/blob/master/quan_matters/test/fusion/inner_product.cpp

I found that compilation was too slow for most real uses ( may be good for rotation matrices etc), but it was fun getting it working

[kwikius]

Added a somewhat more complicated example with vectors
https://github.com/kwikius/units/blob/andy_master/example/trilateration.cpp

This is just a translation of my Trilateration example to mpusz/units
and uses some glue to get mpusz/units working with quan types

requires to include my quan library ( headers only required) to run
Anyway shows a use for unit aware 3d points vectors etc

example output
A = sphere(centre = [0 km, 0 km, 0 km], radius = 8 km)
B = sphere(centre = [10 km, 0 km, 0 km], radius = 5 km)
C = sphere(centre = [7 km, 7 km, 0 km], radius = 7 km)
intersection point = [6.95 km, 1.12143 km, 3.79999 km]

[jml1795]

Was taking a look through both the proposal and some of the discussion. I really like the idea of standard support for a units library and have tried some of the other referenced libraries in both toy and production code in the past.

Regarding LA support for units, one issue could be matrices and vectors that are inherently of mixed dimension. Something like a covariance matrix where you have a position x position block, velocity x velocity block, and pos x vel blocks likely end up causing issues. Not sure there is a work around, but worth giving it a thought.

[hatcat]

The linear algebra proposal, https://wg21.link/P1385, allows for customisation of all operators so you can provide your own algorithms. When you declare your particular vector or matrix type, the operator traits object is part of the spelling.

[mpusz]

Last week in Prague we made Untis library to work with P1385 :-) There are some changes still needed to be done but soon I will push the LA integration and some examples.

Regarding the question about storing different types in one vector/matrix, I am not aware of any LA library that does this (but it does not mean it is a bad idea).

[jml1795]

Things sound promising for the numerics proposals. I have followed 1385 from the sidelines for a while and will read the traits section a bit closer.

Regarding the question about storing different types in one vector/matrix, I am not aware of any LA library that does this

Me neither. Not sure it makes sense to explore deeply or if there is a better place for the question (is this a units question, a LA question, a numerics question?), but the problem is frequently encountered in many fields, particularly in tracking and perception systems.

A simple example could be a state vector from a Kalman filter with position and velocity state elements, possibly units::si::length<units::si::metre> and units::si::velocity<units::si::metre_per_second>. What would the syntax for declaring such a thing be, if possible? I've thought about it in the past as something like a (hopefully contiguous) tuple with the common LA operations defined on it.

[dwith-ts]

Not sure if that was clear from my earlier comment: I have written a comprehensive library that can handle vectors and matrices with mixed units. The library is used in production (unfortunately not open source at the moment).

It is also able to completely handle the Kalman filter for tracking and perception systems example @jml1795 mentioned (that was actually the use-case I developed it for). A state can e.g. contain x- and y position, velocity, acceleration and orientation, the same holds for a covariance matrix.

Design-wise I decided to make both the underlying units library as well as the underlying linalg library exchangeable.
It is possible to get the simple case of a vector / matrix with uniform units to work by customizing the value type of a linalg library (as @mpusz @hatcat and showed). If the general heterogeneous units case should be solved, I found it better to seperate the concerns of

1. providing a units library for scalars
2. providing an efficient linalg library
3. providing a library that annotates vector and matrix rows and columns with units (this is sufficient as all unit-correct matrices are an outer product of a list of row units and column units, see the pretty old book http://www.georgehart.com/research/multanal.html and the short paper www.georgehart.com/research/tdm.ps); this library then uses the libraries from 1. and 2. for its implementation

which also makes each of the parts exchangeable.

From that implementation, I have a few wishes / requests for the underlying units / linear algebra library:

• The units library should be able to handle all scalar types, floating point as well as integral types

• Plane angle (Radian, degree) support would be great (although implementing this correctly is quite challenging and not as straigthforward as one would think)

• The linalg library should have traits that specify whether a type (matrix expression or proper matrix with storage) should be stored by value or by reference in an expression.

• it should clearly state when it returns temporary vectors / matrices instead of matrix expressions (or not do that at all if feasible)

[kwikius]

Here I implemented a matrix using physical quantities, in this case to transform 3d points in length units, but could be of any quantity type
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp

Note that the type of result of a matrix multiplication needs to be deduced. Deducing the types of the elements of rthe resulting matrix takes time,ok for 4x4 matrices but tends to take a long time and memory in the compiler as the size of the matrix is increased
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L433

The basic_matrix class basically just wraps a tuple, that holds the data. The matrix algorithms work on matrix concept rather than type,
https://github.com/kwikius/quan-trunk/blob/master/quan/concepts/fusion/matrix.hpp
so for example a multiplication returns a multiplies view, storing refs to the inputs( which may also be views of course
https://github.com/kwikius/quan-trunk/blob/master/quan/fusion/make_multiplies_view.hpp

https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L117

The positions coordinates is first converted from a 3d vector of length into into a homogeneous 4 row x 1column matrix. The extra point is required for compatibility with 4 x 4 transform matrices, to do perspective transformation etc
https://github.com/kwikius/quan-trunk/blob/master/quan/fusion/make_column_matrix.hpp#L25

functions returning matrices for the common functions are shown ( the Q type is either a physical quantity or a numeric type)
// translation matrix
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L131

// rotation around x
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L160
//rotation around y
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L189
//rotation around z
https://github.com/kwikius/Trilateration/blob/master/trilateration_transform_matrix.cpp#L218

[dwith-ts]

Note that the type of result of a matrix multiplication needs to be deduced. Deducing the types of the elements of rthe resulting matrix takes time...
This deduction using all matrix entry types can be simplified by only using each matrices row + column types, see also the documents I linked.

When calculating A*b, the inner units (A’s column units and b’s row units) have to cancel which e.g. happens if A is a Jacobian (a transformation is the Jacobian of the linear transform between 2 vector spaces and the units in a Jacobian are given by the row unit divided by the column unit) and b is a vector.

These rules are similar to the classical matrix multiplication rules where the inner dimensions have to match. In the cases of Unit-tagged matrixes, the inner units have to be the inverse of one another in order to cancel

[dwith-ts]

The actual calculations can all be done by an underlying matrix library, the responsibility of such unit-safe matrix classes is mainly to promote the types during calculations, sprinkling the syntactic sugar on top.

[kwikius]

@dwith-ts The links on dimensioned matrices are great. I have never seen any literature about using physical quantities in matrices before. If compile time can be reduced dramatically for larger matrices, then that makes the use of such matrices much more practical. Thanks for the information!

[mpusz]

Done: https://mpusz.github.io/units/use_cases/linear_algebra.html :-)

[mpusz]

Please note the issues for LA: BobSteagall/wg21#36, BobSteagall/wg21#37, BobSteagall/wg21#38, BobSteagall/wg21#39.

[dwith-ts]

For your information: I just submitted a presentation about “Physical Units for Matrices” to CppCon2021. If this will be accepted for the conference, I am going to explain in detail how I did the implementation.
As I already explained, the big advantage is that it also works for the case of non-uniform units (i. e. each element can have a different unit)

[JohelEGP]

Looking forward to it. I believe std::bit_cast might be able to do away with concerns about a matrix of quantities/quantity of matrices and their usability with existing linear algebra libraries that may require something like int[3] to take advantage of built-in CPU instructions.

[dwith-ts]

Std::bit_cast could also be of help there I guess, but it seems that this returns by value, i.e. it makes a copy? This would produce additional run-time overhead then. Or do you know a way how the copy can be avoided?

[JohelEGP]

std::bit_cast alone can't help that. Specially when it come to bigger matrices, like 3x3 or int[9].

With platform-specific knowledge, you could turn the call to std::bit_cast into a reinterpret_cast when not std::is_constant_evaluated(). That's a minefield that requires proper exploration, though.

[kwikius]

@dwith-ts. This sounds great. I did implement matrices with arbitrary types which work with quantities . For example
https://github.com/kwikius/quan-trunk/blob/master/quan_matters/examples/fusion/kalman1.cpp#L69
However my implementation was quite "dumb".

The main problem I found was that compile times grew very fast with larger matrices. What are the largest matrices you have used? If you have got an efficient implementation for larger matrices I would love to try it out!

(The matrix '''x''' above actually uses 3D vectors as elements so the result of m1 * x is really a 3 x 3 matrix) I found that anything over a (say) 5 x 5 matrix was impractical in my implementation

FWIW A general ( un unit related) advantage of using a matrix with individually typed elements is that it is possible to optimise out calculations where the result is known at compile time. The above matrix '''m1''' elements are actually so called static values.
These have the property that by design no runtime calculation is performed where the result is known at compile time. For example a multiplication of a runtime value by a static_value representing 0 will return a ( suitably typed) static_value representing 0. Further they only use the minimal one byte for storage ( The type information for arbitrary types includng physical quantities is carried along in the type. The value is held as a compile time rational. In any calculation where the result is not known at compile time or can't be held exactly, the static value is evaluated to its normal runtime counterpart. The type of these elements in the result matrix is of course calculated at compile time as part of the signature of the matrix multiply function.

[dwith-ts]

The biggest matrices I tried so far are 9x9. However, it should work with matrices of any size as there is no need to actually model the unit of each entry and thus having M*N template arguments to a matrix.
All useful multipliable matrices of size MxN have a dimensional freedom of M+N because their dimensions can be described as an outer product of a vector of row units row_units and a vector of column units column_units:
[Mat_i,j] = [row_units_i] * [col_units_j]
(This can even be reduced to M+N-1 as the first entry of row_units can be restricted to be dimensionless while still being able to represent all possible matrix units)
Furthermore, my solution does not try to actually put units inside the matrix, this will IMHO break down pretty quickly as the compiler dies of template overload ;-) Instead, the user-interface of my matrix class (and its type promotions) take care of annotating / returning everything as an appropriate physical quantity.
@kwikius your idea of optimizing out certain operations in such matrices (like multiplication with 0 or 1) is really nice, this is also something that I already tried out and it worked.

[kwikius]

OK. So (something like) the units of elements of an m * n matrix M can be worked out from a row matrix rM and a column matrix cM, because all dimensioned matrices have a specific form as cited in the refs above.

• in pseudoC++ The unit of an element of M[r,c] can be found from decltype(rM[r] * cM[c]) and most importantly this form reduces the compiletime from O(m*n) to O(n+n) or even slightly better.
• The numerical values are held privately in the matrix as dimensionless types and coerced to the correct type when required?
• The dimensioned matrix framework can be layered on top of existing matrix frameworks

I don't know how far I have it right but no matter. This sounds great and I hope that it gets on the agenda at the conference.

As George Hart hints in his papers, there do seem to be many useful ( and under-explored) properties of dimensioned matrices for example in the areas of type checking and correctnesss verifying and there also seem to be some opportunities for optimisation too.

Are you planning to release any source code? I find it hard to follow abstract maths and I would certainly try and get time to work through or just create some examples.

[dwith-ts]

@kwikius your bullet point list is correct. However, if you really want to implement something it IMHO makes sense to wait for my presentation, there are a few other tricks I don't want to spoil just yet (apologies for the secretiveness).

Quoting from my submission overview:
...Furthermore, the approach provides an even richer (and more type-safe) annotation for linear algebra types and their elements than only physical units: it enables annotations that resemble quantity kinds. This means that it can also distinguish between X- and Y-position in the same coordinate frame and between X-position in coordinate frame A and X-position in coordinate frame B. In addition to the element annotations, we also get stronger matrix types which encode the content, e.g. a jacobian matrix, covariance matrix, position and displacement vector types. I will show how these annotations propagate through linear algebra operations, how they determine the subset of valid operations on each type and of course and most importantly, how this can be implemented efficiently in C++. This will enable a whole new level of type-safety in C++ linear algebra applications that allows to catch the majority of mathematically inappropriate operations on vectors and matrices, whereas normal linear algebra libraries like e.g. Eigen, blaze or P1385 can only catch problems related to incompatible vector and matrix sizes. ...

@kwikius : I am trying to get permission to release source code, but I do not want to promise anything here.

[mpusz]

This sounds great! I always hoped that someone will pick up this subject and provide a standalone library on top of mp-units and linear algebra proposals that will merge them together in an efficient way.

I hope to see you at CppCon. We definitely have to meet over lunch there and talk about it 🍽️ 🍻

[dwith-ts]

I hope to see you at CppCon. We definitely have to meet over lunch there and talk about it 🍽️ 🍻

Perfect, let’s do that! Looking forward to discussing this with you! 🍻

[dwith-ts]

FYI: I am going to present most of my submitted CppCon talk at MUC++ next week https://www.meetup.com/de-DE/MUCplusplus/events/279334186/ , you can join after signing up. Hope to see a few of you there!

[JoostHouben]

@mpusz this issue was closed with a link to https://mpusz.github.io/units/use_cases/linear_algebra.html, but that link is now dead. I could not find any current documentation on linear algebra support with mp-units, but perhaps I overlooked it? What is the status of linear algebra support in mp-units currently? Is there any inter-operability with Eigen?

[mpusz]

Hi @JoostHouben! The problem with the linear algebra support is that none of the libraries on the C++ market does it right. And the only library outside of C++ that I am aware of is Pint. By this, I mean that none of the physical units libraries allow using a linear algebra type as a representation type and doing correct vector math on it (dot/cross product, etc) to obtain proper quantities.

In 2.0, we are refactoring the previous support, and it is still WIP. You can find some documentation already here: https://mpusz.github.io/mp-units/2.2/users_guide/framework_basics/character_of_a_quantity. You can also find some code examples here: https://github.com/mpusz/mp-units/blob/master/test/unit_test/runtime/linear_algebra_test.cpp.

There is also a PR #493 started, but I did not have time to work on it for a while.

All of those are subject to change. Early feedback is really welcomed as well.