Linear algebra phone meeting 03 Apr 2019 Present: - Bob Steagall - Guy Davidson - Michael Wong - Steven Varga (Toronto) - Mark Hoemmen - Cem Bassoy - Graham Lopez (ORNL) - Christian Trott - Jayesh Badwaik MW: SG14 will not meet next Wednesday, due to the ACCU meeting in Bristol. Will move SG14 meeting to the followingx week. BS: Progress report on P1385. In our public GitHub repo, we have a new branch "la_r2" containing work in progress. I've been implementing KONA 2019 suggestions, documented at end of P1385R1 (at wg21.link). I've removed column/row_vector and gone to a single vector type. I took a detour and started working on something called "number traits" -- an initial attempt to differentiate between semirings, rings, and fields. Matrix elements would need to be fields. Not sure if that's the right thing yet, but I wanted to put the machinery in place. I've also taken dynamically resizable etc. engine classes, and split them into separate classes based on whether they support matrices or vectors. Not sure if that's the right design, but I did it as an experiment. I may instead partially specialize on the number of dimensions. BS: After this, I will look at integrating mdspan. Finally, I start thinking about executors and execution contexts. It will likely be several weeks before I get to these two tasks. MH: Regarding semirings, rings, and fields: Different algorithms have different requirements on the matrix element type. For example, algorithms that look like tensor products (e.g., matrix-matrix multiply, dot products) only need plus and times. Factorizations may need division, but in some cases could be reformulated to work without it -- e.g., for matrices whose elements are polynomials. Would it make sense to let each algorithm constrain the matrix element type, or do you want to constrain the matrix element type for all algorithms? BS: I was thinking about that with the traits refactor. Naively, it suffices for the matrix element type to be a ring. But what about the determinant? That needs division. For now, the test is that the matrix element type must be a field. MH: Would you prefer to express this requirement as a Mandates (static_assert) or as a Constrains (SFINAE, that is, does not participate in overload resolution if requirement not met)? BS: I prefer mandates (static_assert). MH: I would have a preference for a per-algorithm requirement, and to express the requirement as Constrains (SFINAE). However, our use cases for ring or semiring matrix element types are mainly for graph algorithms, and would thus use sparse matrices. MW: Is this work on semirings, rings, and fields for the first layer? BS: If by "first layer" you mean the foundation layer, then yes. It's something new I thought of a few days ago, and added just to start the conversation. Not sure if these traits should be part of linear algebra, but if useful there and applicable elsewhere, then it may make sense to pull them out into their own paper, for Numerics. MW: OK. We should check thoroughly how much we need that for first layer. MH: ... field theory ... [I think I said something to the effect of "this is interesting and potentially useful, but we need use cases for more exotic matrix element types, and we also risk accusations of overcomplicating the design to support them." In addition, interfaces of libraries like the BLAS and LAPACK were not designed for this much generality, so we may need to look elsewhere to make sure that the interface itself does not overconstrain us -- especially as we start to include matrix factorizations.] MW: I see what you mean. GD: When I started working on this, I wondered how far to draw the line [with respect to how general the matrix element type should be]. I do use matrices of polynomials. BS: At least from an interface and type design perspective, having a set of elements that's well behaved and fulfills a fixed set of requirements ... [is a good idea]. We could have an element type that represents a polynomial, that has well-formed operations ... MH: This issue of the matrix element type matters more for matrix factorizations, than it does for algorithms that look like tensor products [that only use plus and times]. BS: Any more comments on P1385? MW: Will you submit a revision for the pre-Cologne mailing? BS: Yes, we would submit an R2 paper for the pre-Cologne mailing. MW: Someone from PNNL, working w/ Andrew Lumsdaine, wants to give a talk at the SG14 linear algebra meeting in Cologne. I would book a face-to-face meeting at Cologne. GD: I will be there. BS: I will try to be there. === BS: Are people planning succession papers [to P1385]? GD: I've started trying to work out what a geometry succession paper would look like. I started with Boost Geometry and with our games, to see what minimum set would be, so that others could write further succession papers. MH: Geometry is numerically subtle. GD: It is. MH: For example, correct evaluation of Boolean predicates (like "am I inside or outside of this circle?" or "on which side of a line am I?") may call for extended-precision arithmetic. GD: Games usually express this with an "epsilon" tolerance, that lets them trade off between speed and accuracy. MH: Finite-element or other discretizations of partial differential equations in space favor accuracy, in order to avoid situations like inverted elements that may break the discretization. BS: This reminds me that I'm currently using a geometry library for earth surface stuff, open source, originally Java, then ported to C++, not well. But it's mandated that we must use it. MH: Regarding succession papers: We're working on identifying "primitives" for a linear algebra library -- low-level things in support of an interface like P1385, that vendors would most want to implement for good performance and/or accuracy. I am working on a paper that shows how a C++ linear algebra interface might develop organically, starting with a BLAS interface and developing to higher levels of abstraction. The paper is not ready yet, but I could summarize it now, if you like. BS: Please discuss it at the next call, once it is complete. MH: Regarding the set of supported matrix element types: I'm worried that we're taking as models interfaces that only support floating-point types. Would we need different interfaces to support fixed-point types? [For example, would they need to pass along scaling factors?] We should consult experts on fixed-point arithmetic. GD: I'll try to get John Mcfarlane on the next call. GL: Davis Herring at LANL has interest in fixed point and also in linear algebra. BS: Davis was in the LEWGI review of P1385 at Kona. GD: He took excellent minutes. GL: He was representing one of the fixed-point papers as well. === BS: Misc. business -- other topics CB: Google Summer of Code started. I am contributor to muBLAS and added tensor library to that. Next GSC topic: I want an interface to that, similar to your paper. I asked Jayesh if he wanted to participate and he said yes. I would like to extend this ... to tensors. BS: I'm interested in that ... based on the work we do here, some subset of us would like to put together a tensor paper proposal, once we've learned lessons from linear algebra. CB: We're trying to use Boost expression template evaluator -- can do Eigen Tensor - like things -- specify at end of expression whether you want to evaluate in parallel, on device, etc. Trying to do that in this Google Summer of Code project. We have several good students. MW: Would the tensor paper you're suggesting also be of interest to SG19 (machine learning)? CB: Yes. In Boost, we are heading twoards tensor decomposition algorithms, like higher-order SVD. In machine learning, you may have other operations but similar. That would be a good thing to do. MW: Christian Trott is here with me; add him to roll call. SV: As I mentioned [at the beginning of this call, I work on a ... tensor solution for the HDF group, for linear algebra library, ... idea is general framework for C++ ... current status is full linear algebra support ... structured ... compiler-assisted reflection ... how does this fit into the material that I see on my screen? and how related to machine learning working groups? ... BS: I don't know yet. I'm looking at your link. SV: Is it OK for me to attend these meetings to get a better understanding of the trajectory towards which you are working? GD: That would be great. These are public meetings. You don't need permission to join us. BS: You're welcome any time. MW: SV, since you live in Toronto, you can apply to be a Canadian expert -- just send me an e-mail with your resume. With that, you would have the power to vote officially. Every country has its own delegation -- SCC -- I am the chair of all the programming languages. SV: One more question: I just received this material from you, BS, are these functions functional? If they are working, then I would incorporate them. BS: What we're doing with P1385 is to specify an interface that would become part of the Standard. Library implementers would provide their own implementations that express the interface that we propose. [Jayesh Badwaik joins] BS: [continues] Our work so far, up to this point and probably for the next few weeks yet, is on the interface, and perhaps a reference implementation that vendors could use or modify. Most of our efforts now are on designing the interface. Branch la_r2 that I sent you -- la_r1 branch actually compiles and executes, but with stub execution. Once I get r2 branch working, the same thing [stubs]. Beyond that, we would actually put real code that does arithmetic there, into the arithmetic traits where the work would actually be done. Also, this interface is deliberately limited to what you're calling rank 1 and rank 2 objects -- matrices and vectors. BS: While working on traits for semiring, ring, and field -- several people asked me the same Q about element promotion traits -- they asked why we couldn't just use common_type. At the time, I didn't have a good answer. common_type is one of those things in namespace std you're allowed to specialize. Half-size floats: It might be that you want to write traits for operations, where both operands have type half-size float, but you want result type to be full-size float. It could be that common_type of half-size and half-size is half size. MH: common_type of complex and double is complex; surprise! GD: Wow. BS: I'm writing a paper to fix that aspect of common_type; Walter Brown discouraged me from doing so, though. [The implication to the scribe was that WB thought it would be impossible to get this change through, not necessarily that WB thought it was a bad idea. I'm not sure exactly what BS meant, though.] GD: WB also discouraged me. MH: Would the expression C = A + B for half-precision A and B and full-precision C trigger half,half,full promotion, with correct arithmetic that accurately computes a full-precision result? or would it just truncate half+half to half precision and cast to full? BS: The intent is that you could support the former, yes. MW: You have two meetings before the mailing deadline -- May and June. CB: Would you like to discuss about the view stuff? BS: Yes. CB: [shows an example of views] In a Householder QR decomposition, in Matlab you would use a colon operator. Here, I show "span(j,m)". The idea is to create a matrix view. No copying, it would just be a view. Do we want to support something like this in the matrix library? BS: It's not part of the interface now, but in KONA feedback, there's a desire to have view types of columns, rows, and submatrices, that are part of the library. We haven't defined them yet, or decided whether they would use mdspan or something else. Christian Trott: Another argument in discussions with vendors why we should use something like mdspan -- I talked with vendors about accessors that express things like streaming stores and loads -- using special capabilities in the memory subsystem -- mdspan has all the extension points needed in order to express all these things. Vendors are very interested in that. On newer systems, memory system hardware is more complex and can express more than just plain load and store. CB: Of course mdspan supports a lot of things, but if you want to keep it a bit simpler ... [shows screen] -- operator() that takes "spans." MH: We have a function called "subspan" that does this. Can be hard to have overloaded operator... would be good to use mdspan here since it has gone through more design review and more C++ subcommittees, but we don't mean to stomp on people. GD: No problem from us -- we view mdspan as a crucial implementation detail, potentially as a view engine, so please don't worry. JB: Sorry I haven't had so much time to work on this; I'm working on my thesis. [group]: Don't worry! === BS: Next linear algebra call: Wed 01 May. MW: Next SG14 call: 17 Apr, 2-4 PM Eastern.