# C++ linear algebra phone meeting Date: 2019-06-06 15:00-05:00 Leader: Bob Steagall Scribes: Mark Hoemmen, Bob Steagall ## Participants - Jayesh Badwaik (Jayesh) - Charles Bay - Matthew Butler - Guy Davidson (Guy) - Ronen Friedman - Mark Hoemmen (mfh) - Klaus Iglberger (Klaus) - Graham Lopez (Graham) - Maikel Nadolski (Maikel) - Rene Rivera - Brett Searles (Brett) - Bob Steagall (Bob) - Michael Wong (Michael) Bob: Update on P1385. As discussed in May, Guy and I made changes to improve customizability. We presented that at C++Now in Aspen. I got good feedback. Here is a summary: One-vector approach (that is, not having row and column vectors be separate types) turned out to be a controversial design decision. Audience favored using the std::math namespace. Regarding the traits-based approach: I got feedback that it seemed like a "1990's approach" to writing an interface. The commenter didn't suggest a way to implement this functionality without traits. However, I talked with Zach Laine who did outline an approach. The idea is to remove the second template parameter (operation traits) from the matrix and vector class types. Positive feedback on customizability -- people liked arithmetic traits -- some discussion of predicate traits, e.g., number traits for ring vs. field, etc. We were asked why we didn't just use concepts. The answer is that we haven't had time to figure out how to use concepts to do that yet. I had a conversation w/ Lisa Lippincott. She liked the analysis that preceded the actual design. I had a conversation w/ David Hollman (Sandia). David mentioned the papers Sandia is developing. Guy and I are currently working on a revision of P1385 for pre-Cologne submission. --- Bob: Klaus has done a nice writeup of feedback for P1385. Klaus: We [the proposal?] do not deal w/ complex expressions. It is possible w/ infix [arithmetic operator] notation to combine arbitrary number of vectors and matrices [in an expression]. How do we want to deal with this? e.g., A+B+C is easy. but (A+B)*C ? This makes us decide how to deal with the addition. There are 3 different approaches, each of which has been taken by different libraries: 1. Do nothing / forbid this expression from compiling (Eigen) 2. Treat A+B as an expression (template), so that the matrix-matrix multiply evaluates A(i,j)+B(i,j) every time it needs that entry of the expression. 3. Create a temporary matrix for the matrix sum A+B. We must decide whether or not this expression is allowed. We can't let whether or not it works be implementation defined. Bob: I don't have a direct answer for that in terms of P1385. Guy and I hesitate to prescribe how to evaluate expressions; I would rather call this "implementation defined." Maikel: "You want to specify order of evaluation more strictly?" Klaus: I don't think you should define the underlying technology, but you specify whether it's possible. Regarding Maikel's question, it's not about evaluation order, but whether to allow certain expressions. For example, (A+B)*C is a valid compound expression, but should we allow it? How it works is not what we should specify, but whether it works is. Klaus: Also, should taking a submatrix of an expression, like "submatrix(A*B)," be allowed? or should "submatrix" only work on a concrete matrix type [not an expression]? mdspan may not be enough to support this use case if we need taking a submatrix of an expression to work. mfh: It's possible to do some expression templates w/ mdspan through its Accessor policy; the question is whether you want to do expression templates. Klaus: That's right. Klaus: Last time, everybody wanted begin and end (iterators) for vectors. However, they did not argue for iterators for matrices, because this is more difficult. But you'll need them once you introduce sparse matrices. Should we deal with this now or wait until add sparse matrices? Once we have iterators on matrices, should we introduce Standard-compliant begin and end member functions that take no arguments, or have them take an index argument and thus allow traversal of a single row or column? The former (zero arguments) has subtle implementation disadvantages: more memory, larger iterator, and slower. The latter would not be STL compliant, but it might work better for matrices in general. mfh: mdspan only allows indexing; it does not offer iterators. This is a deliberate choice. Klaus: But sparse matrices. mfh: The issue with sparse matrices is that users are less likely to want to construct many different kinds of arithmetic expressions with them. However, I can suggest prior work -- in C++, even -- on using relational algebra to define something like iterators over different sparse matrix data structures, for the automatic construction of efficient computational kernels like sparse matrix-vector multiply. See "The Bernoulli Generic Matrix Library" and other publications by Nikolay Nateev, Keshav Pingali, and Paul Stodghill. Brett: Could we modify Ranges for spanning matrices? I know it "uses iterators" yet it has better capabilities. mfh: Ranges are 1-D, just like iterators. Klaus: Regarding the examples of static and dynamic matrix type promotion in Section 6.6 of P1385: The entire text focuses on multiplication. For addition, dynamic + static could perfectly well result in a static matrix. Guy has talked about this at Meeting C++, but it is just not described yet in the paper. Feedback: Be more specific about how result types are evaluated, and give more examples in this section. Guy: You're correct, we did talk about this at Meeting C++, but only multiplication. Addition is different because dimensions of output are same as those of input. I haven't given thought to the other operations other than multi and their effects on the return types. I will think about that more; thank you. Klaus: Might need to introduce operations that don't make sense in generic case, like e.g., Kronecker product. If we have an operation for matrices, we should also think about whether it makes sense for vectors. Klaus: I have 6 pages of findings that I can send around. Guy: Thank you, Klaus. Klaus: Any more questions? Bob: Klaus, you expressed an opinion in favor of row and column vector being separate types. I'd like to hear your thoughts. Klaus: I think it's valuable to have both, in my experience. I'm happy to distinguish between the two. I deal with more complex stuff, like treating a vector like it's a matrix -- for that, it's helpful to know whether it's a row or column vector. Also, sometimes we want to treat a vector as a matrix, or a matrix as a vector. I like the fact about math that it's exact, but here we're being vague -- we're interpreting something because we feel that it's what users want. Graham: This opens a can of worms... Bob: Please go ahead. GL: 1. We want to have an interface... it's impossible to cover every edge case, so we want to cover the most common but allow special cases. For example, if you want to treat a vector as a single-column matrix, you already have the machinery to do that. We can provide helpers to do that cheaply. But, it's much more often that you want to treat a vector as a vector, not in the linear algebra textbook sense (which is the only place where you differentiate between row and column vectors). 2. Treating a vector as a matrix is type punning, which is generally not accepted for interface design. I'd like to call types what they are. mfh: Our job is to create low-level machinery, not to please users w/ terms... [goes on to explain that decisions like whether row and column vectors should have different types, or whether vectors get automatically promoted to matrices or vice versa as needed, are high-level decisions on which users disagree.] Bob: The point of 1385 is to create basic facilities and interfaces for users, and let us gain experience w/ that, and to add more sophisticated things later. I would like to punt on the idea of how to manage expressions. We don't have to handle that in this paper. Could leave for follow-on papers. mfh: Precedent for that in valarray. [valarray permits but does not require expression templates. There is a bit of history there. Daveed Vandevoorde attempted to add expression templates to the valarray proposal, which had been abandoned by its original authors, but it was too late in the standardization process and the Committee instead accepted the current language. Some implementations of valarray do actually use expression templates and have other optimizations.] Bob: We did start out with a small set of requirements we were trying to fulfill. Our paper tries to fulfilll those requirements. We can ask whether they are the right set of requirements. Graham: I agree with trying to scope it down to provide the low-level operations. That's gives us a smaller goal and less to argue about. I would like to avoid trying to make things nicer for the user prematurely. Nevertheless, I would be open to arguments about the two-vector system, if arguments could be made for it independent from usability. Are there performance implications? Could you not express some operations in the one-vector system, that you could in the two-vector system? Klaus: That would only come up later (more operations, not defined in the paper). I think there will be second, third, and fourth steps. By saying there is only one kind of vector, do we restrict ourselves in the future? Bob: Going back to a conversation Graham and I had in Kona: my original thinking for row and column vectors was to distinguish bewteen inner and outer products when using operator*. Maikel: My experience from using Blaze or Eigen is that differentiating between row and column vectors caught a lot of bugs at compile time. Jayesh: I agree, esp. in generic physics code. Graham: Outer and inner product are not the only possible vector products. Row and column vector thus won't get you all the way. Bob: Any more comments? [none] Bob: Klaus, thank you. Turn over to Jayesh. --- Jayesh: [shares screen] I was wondering how to deal w/ problems like what Klaus raised today, w/ submatrix(A*B). Also, as we put stuff into the Standard Library, we want it to be generic and customizable -- this puts a lot of pressure on the library. Is there a way to let users extend the Standard library, or to use their own library w/ Standard Library types? e.g., std::sort or boost::sort w/ custom containers. Problem is not the named functions themselves, but with infix operators like "a+b". Is this an expression template or an owning type? What if a and b belong to different libraries? Here is an example. Sorry for "engine" overloaded term w/ P1385. Please don't confuse w/ how I'm using the term. Also, this could be written in terms of concepts, but I'm not so familiar w/ concepts so I'm writing using traits. [presents design; the idea is to distinguish between owning types, consuming owning types, and expressions (that don't own). Define operations on engine types. Objects like vector can change engine. operator+ defined on T&& and U&& but constrained so they use a particular engine. Example shows how to add two vectors, one that uses a Standard Library type and the other that uses a custom type.] Jayesh: I have a proof of concept that the two mechanisms work, and am currently trying to combine them. ... if we place similar requirements on our vector types, e.g., that they should have an Engine template parameter, export types, and have the three change_engine functions, then we can get interoperability between libraries. [Around this time, Guy starts having internet trouble and dropping out.] Jayesh: ... and we do not have to restrict how the execution of the addition is being done. e.g., if you want to do a different addition for sparse matrices... Bob: Feedback? mfh: Sparse matrices are a different world. The kinds of operations users want to do with sparse matrices is often very different than with dense matrices. Thus, you don't have to feel like you need to solve that problem. Jayesh: Point is to let the Standard Library move slowly.... Bob: Other comments or questions for Jayesh? mfh: Just wanted to say that this is good stuff! It looks like this could help decouple use of arithmetic operators from linear algebra. Michael: Jayesh, will you submit this to Cologne? Jayesh: Not sure; depends if this paper is ready. It is very preliminary. However, I will be there in Cologne. Michael: I will try to schedule a [face-to-face] SG14 meeting in Cologne for Friday morning. Bob: Jayesh, this is interesting work. Even if you can't finish the paper in time for the deadline, please keep us appraised. --- [At this point, mfh talks about the 3-4 papers on linear algebra that Sandia intends to submit.] [Rene Rivera and Charles Bay enter the meeting.] --- Bob taking notes now ... mfh: several papers under development. * md_array proposal, like span, but an owning container like std::vector. Will let small data structures be passed around by value. * another papers that describe lowest-level things that could possibly work. Started with idea of building upon BLAS; an experience paper describing how to wrap BLAS in C++; another aspect is how to simplify how arguments are managed * another paper is how to support batched linear algebra operations; * papers still under review, not quite ready for publication * does not compete with 1385, but provides a set of "atoms" upon which 1385 might be built. * may be able to present all 3 at Cologne * not attempting to provide expression support * could do simd inside Michael: To the group: Is it appropriate to shoot for a TS? Bob: Not sure; whichever approach gets linear algebra into '23 fastest is what I want. mfh: There seems to be discouragement against using TS's lately, don't think it's the right way to go. Michael: Not so much for Library, but agreed, a TS may be too much. Michael: Do you intend for your proposal to be a TS? Bob: no mfh: no --- mfh taking notes now --- [We interpreted general feedback from previous meetings as discouraging a TS.] Michael: SG14 monthly call next Wednesday; it will be the last before 17 June pre-Cologne mailing deadline. Bob: 2019-07-03 15:00-05:00 will be next call.