Different types for matrix columns

[mpusz]

Could you please consider specifying different representation types for matrix columns? It would address a lot of engineering needs. For example, together with the physical units library, we could implement matrices as described here: https://www.kalmanfilter.net/default.aspx.

[mhoemmen]

Commenting because I'm not quite sure what @mpusz means and would like to learn more.

[mpusz]

I would like to create a matrix like below where L, V, A are values of distinct types:

| L V A |
| L V A |
| L V A |

For example, L could be length, V a velocity, A an acceleration and together they could be parameters of a quadratic equation:

x = x0 + v0 Δt + 12 a Δt^2

I got a lot of requests like the above recently to be supported by my physical units library (https://github.com/mpusz/units), but it actually is a requirement for the LA library as my library will work with such a library just fine.

[mhoemmen]

@mpusz thanks for explaining! Would you object to the library implementing this by storing type information separately from a matrix of unitless quantities? For instance, one could treat the matrix in your example as a 3x3 matrix of unitless quantities, times a diagonal matrix diag(1.0 * L, 1.0 * V, 1.0 * A) (sorry for the Matlab-ish notation; could really use math LaTeX here). That would be helpful for keeping the units around while still being able to call third-party libraries like the BLAS or LAPACK.

[mpusz]

I think it would work but we still need to be able to create this diagonal matrix/vector storing 3 different types, right? If we do it, I think that it could be easily extended to regular matrices?

[mhoemmen]

I don't want to hijack Bob's issue tracking system ; - ) . The idea I had in mind is that there's a performance and build size/time cost to instantiating operations like matrix-matrix multiply with arbitrary types, so I'd rather try to "algebra out" the units and just feed unitless quantities into a low-level "math kernels" library.

What I'm trying to say is that I don't object to your idea; I just might choose to implement it by storing the quantities separately from the units (sorry if I'm getting the terminology wrong). On the other hand, I've spent a nontrivial part of the last 10 years of my career reducing build sizes and times by controlling instantiations, so I'm probably a bit biased against pushing the type information into the lowest-level representations.

[BobSteagall]

Hi Mat and Mark. I'll have to think about this. At first glance Mark's diagonal matrix idea seems plausible.

However, one barrier is the underlying assumption in the design that all elements of a matrix/vector are of the same type. No thought has been given to rows/columns of different types.

It may be achievable with custom engines and custom arithmetic traits, but would require some dedicated time to thinking about it.

[mhoemmen]

@BobSteagall Hi Bob! :-D My feeling is that this is something you would want to lift above the level of "element type," because the kinds of things people like to do with matrices would constrain the units of the elements of a matrix in ways that might be hard to check at compile time. For example, would I need to "solve for the units" when solving the linear system Ax=b ? Do I need to do that at compile time, if my interface looks like x = solve(A,b)? I haven't thought about what that would look like, but it's unprecedented (as far as I know) in our world of large sparse matrices to encode units in the types of the matrix elements at the lowest level of representation.

[BobSteagall]

Mark, I agree. One can imagine a kind of "shadow matrix" that carries the extra type information (beyond unitless element type, that is).

[mpusz]

it's unprecedented (as far as I know) in our world of large sparse matrices to encode units in the types of the matrix elements

Yeah, I agree, this is a current state. However, it is like that only because until now we did not have a good modern C++ LA library and standardized physical units library. It might change with C++23 (or C++26) and this is why the first requirements for such use cases are raised to me by the users and I pass them to you as they require a modern C++ LA library ;-)

[mhoemmen]

Hold on, let's not talk about "modern C++ LA libraries" quite just yet. Doing this efficiently will take some thought; it's not just a "multidimensional tuple."

Our users use our C++ linear algebra library (not P1673; a much more complicated one) to solve linear systems with dimensions too big to fit in 32-bit integers. It has millions of lines of code, and implements a variety of algorithms, from algebraic multigrid to complete and incomplete sparse factorizations to Krylov subspace methods. Our users care about units, but they can and need to control them at run time. They will not accept a library that needs to be reinstantiated for every possible permutation of units. This means that some level of the implementation will need to strip away units-related compile-time type information.

I welcome thoughts about how to do that correctly, but I don't like the implication that "modern" necessarily means "done at compile time" or "pushed down to the lowest level of abstraction." Hardware exists that implements linear algebra operations directly as machine instructions, and this idea is becoming more popular with all the cycles that companies spend on machine learning. If we aim for "zero-overhead abstractions," then we need to design libraries that can get matrices, vectors, etc. down into the hardware as bulk homogeneous data. Our C++ libraries need to be able to call optimized low-level functions that vendors have already implemented or will implement.

I want units too, don't worry! : - ) We just need to think about how to do this right.

[mpusz]

Our users care about units, but they can and need to control them at run time. They will not accept a library that needs to be reinstantiated for every possible permutation of units. This means that some level of the implementation will need to strip away units-related compile-time type information.

I never said that your users are wrong ;-) If this is their use case than fine. They should be allowed to do it. However, there are other users having different requirements and use cases. One of them is to create small and simple 3-dimensional vectors and matrices of different quantities.

I do think that both use cases can easily coexist without any performance hit on any of them. In the worst case, a partial class template specialization will be needed which is not a huge problem.

[dwith-ts]

I wholeheartedly agree with what @mhoemmen suggested, keeping the linear algebra in the low-Level and only enriching it with unit information in a higher level (a wrapper class around the linalg types) makes a lot of sense and for more complex cases than uniform units in a vector / matrix it is the only feasible way.

@mpusz in your example with different column types, it helps to ask what an access like mat.column(idx) should return: idx can be a run-time value which means there is no single return value for such a method, it would have to return a variant to be able to store the different possible column types.

Only when a column access would be done with mat.column< Idx >() where Idx is a constexpr integer would it be possible to return the concrete column type depending on the given idx.

Units for a Matrix are always determined by an outer product of row and column units of the matrix (I‘ll Post a link to a paper with proof in the next post) => it is sufficient to store these row and column units (e.g. 8 units for a 4x4 Matrix) which can be stored as template arguments, so no need for an additional Matrix for the units.

I already answered a similar issue in @mpusz Unit repo, I‘ll add that in the next post as reference.

[dwith-ts]

Here is the post from @mpusz unit library Repo:

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

• providing a units library for scalars

• providing an efficient linalg library

• 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.

[dwith-ts]

Here is the post from @mpusz unit library Repo:

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

• providing a units library for scalars
• providing an efficient linalg library
• 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.

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 in my library.
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).

[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 you there!

[dwith-ts]

The talk from MUC++ can be viewed online and there is also a more recent one from Meeting C++ https://m.youtube.com/watch?v=4LmMwhM8ODI

[chiphogg]

Should this issue be closed as "won't fix"?

It feels as though there is consensus on the idea that the units should "sit on top" of the linear algebra, and that the latter should deal with plain, uniform types. Besides the theoretical aspects, I think @dwith-ts has provided ample evidence that this strategy also works very well in practice.

Is there another use case for supporting different types in matrix columns that I've missed?