std::linalg: Linear Algebra Coming to Standard C++ - Mark Hoemmen - CppCon 2023

cppcon.org/
---
std::linalg: Linear Algebra Coming to Standard C++ - Mark Hoemmen - CppCon 2023
github.com/CppCon/CppCon2023
Many fields depend on linear algebra computations, which include matrix-matrix and matrix-vector multiplies, triangular solves, dot products, and norms. It's hard to implement these fast and accurately for all kinds of number types and data layouts. Wouldn't it be nice if C++ had a built-in library for doing that? Wouldn't it be even nicer if this library used C++ idioms instead of what developers have to do now, which is write nonportable, unsafe, verbose code for calling into an optimized Fortran or C library?
The std::linalg library does just that. It uses the new C++23 feature mdspan to represent matrices and vectors. The library builds on the long history and solid theoretical foundation of the BLAS (Basic Linear Algebra Subroutines), a standard C and Fortran interface with many optimized implementations. The C++ Standard Committee is currently reviewing std::linalg for C++26. The library already has two implementations that work with C++17 or newer compilers, and can take advantage of vendor-specific optimizations. Developers will see how std::linalg can make their C++ safer and more concise without sacrificing performance for use cases that existing BLAS libraries already optimize, while opening up new use cases and potential optimizations.
---
Mark Hoemmen
Mark Hoemmen is a C++ software developer with a background in parallel computing and numerical linear algebra. He joined NVIDIA in spring 2022, and works remotely from Albuquerque, New Mexico, USA. He contributes to various open-source C++ libraries, including CUTLASS, a CUDA C++ library implementing high-performance matrix-matrix multiplication (GEMM) and related computations.
Mark finished his PhD on "communication-avoiding" linear algebra algorithms in 2010. After that, he worked for ten years for Sandia National Laboratories, where he did research on communication-avoiding and fault-tolerant (not at the same time, thankfully) algorithms, and contributed to several scientific computing software projects. He then took a position at a private company, Stellar Science, for two years, before moving to NVIDIA.
Mark's preferred programming language is C++. He has been writing it professionally for 23 years, and has been contributing to the C++ Standard (WG21) process since 2017. Mark is main author of the C++ Standard Library proposal P1673, a linear algebra library based on the BLAS (Basic Linear Algebra Subroutines). He is also coauthor of mdspan (which is in C++23), mdarray, and other related proposals. After C++, he feels most comfortable working in Python, C, Fortran, and Matlab. Mark is familiar with several shared- and distributed-memory parallel programming models, and interactions between them (e.g., CUDA and MPI).
---
Videos Filmed & Edited by Bash Films: www.BashFilms.com
KZbin Channel Managed by Digital Medium Ltd: events.digital-medium.co.uk
---
Registration for CppCon: cppcon.org/registration/
#cppcon #cppprogramming #cpp

[@a_s_tarantoga]

Nice to have this inside the STL. What puzzles me is while having the Eigen library around for years that is providing sensible mathematical syntax through operator overloading, why creating an STL extension with clumsy function names as "matrix_product", "scaled" etc. I hope this proposal takes some more rounds with updates of the syntax...

[@niklkelbon3662]

This is not a problem, in fact we need to write all operation in one expression (for optimizing for example) and we should views instead of matrix types, so operator* is strange here. I think best scenario is something like this: using namespace std::matrix_placeholders; constexpr auto operation = ^(A * B) + C; operation.apply(a, b, c); But its just my first idea

[@bernhardmanfredgruber7156]

The proposal is already in wording review, so AFAIU the design phase has completed. It's meant to provide lower level primitives on which e.g. Eigen could be built. Mark explains this partially during the Q&A at 1:10:10. And if you already like Eigen, you can just keep using Eigen.

[@franciscogerardohernandezr4788]

exactly, even though I am not an Armadillo fan, it has a sweet syntax similar to Matlab, easy to represent mathematical expressions with. Blaze (my choice) follows a similar approach.

[@PixelPulse168]

it's because operator overloading could be confusing. For example, M * A, it could be matrix multiplication or scalar multiplication. Explicit definition improves readability, helps debugging. But I understand your point, it makes code more verbose.

[@luciusvorenus5242]

Because this library is based on BLAS, a low-level standard for linear algebra that has been around much longer (somewhere around 80-x) than the high-level Eigen, and is well designed from an optimization point of view. For example, there are interfaces there for many "patterns" that arise when working with matrices, serial/parallel implementations of which are well described and studied, and these patterns would be impossible to implement in the A=B*C syntax. Although Eigen is much more convenient, I can understand the logic of the developers of the standard who try to rely on something more lover-level. In Eigen, by the way, you can use different BLAS implementations (-DEIGEN_USE_BLAS) for matrix operations used on the backend, I assume with C++26 they will add something like -DEIGEN_USE_STD.

[@tux1968]

This speaker did a spectacular job. Very informative.

[@MarkHoemmen]

Thank you for the kind words! : - D

[@drmikethom]

BLAS (linear algebra) operations are available in optimized forms for each architecture. More complex lin algebra algorithms are being implemented in large frameworks (like PETSc, Trilinos, Hypre, SuperLU, etc.) If C++ coding can provide optimized binaries for these solvers, then there is a great reason to pursue these capabilties.

[@perfectionbox]

Great,now we can consolidate all the vector and matrix types every graphics/physics library includes so we don't need to have 20+ versions.

[@yellowajah]

so the solution to the problem of having a difficult ecosystem is to integrate said ecosystem into the core of the language? Sounds like a downsteam hack to an upstream issue...

[@xenap7001]

@@yellowajahSounds like every language.

[@robertolin4568]

Nah. Now you have 21+ versions instead.

[@aniketbisht2823]

Although std::mdspan took a long time to be standardized, its a really nice abstraction. And one beautiful thing about this library/proposal is that it mostly consists of free functions so that the standard library implementers can keep optimizing them as time goes by because they wouldn't be constrained by the ABI (as they are just functions).

[@policyprogrammer]

You can understand the struggle. After all, who, before the C++ standards committee has ever wrestled with trying to create a robust, efficient, and easy-to-use matrix + linear algebra library? C++, as usual has to be a pioneer! Oh, wait, I'm getting a call... it's from the hundreds of such libraries in dozens of languages that apparently work fine and which almost certainly will outperform whatever makes it into the standard.

[@rasitsimsek9400]

In C++ I'm using the Eigen library.

[@literallynull]

based

[@joopie99aa]

What is the expected experience for end users? Will end users call std::linalg directly, or will they still use 'level 1' libraries like Eigen, which will then wrap std::linalg? If the latter, what functionally changes? Eigen already dispatches to BLAS, and BLAS is already standardized, so it seems this just introduces an extra layer with no real gain?

[@peterevans6752]

Dense material! I like the approach.

[@sucim]

This is a great interface, resonates a lot in how I implemented the dispatch in my own library.

[@bernhardmanfredgruber7156]

I really enjoyed the history section on how the BLAS came to be!

[@MarkHoemmen]

Thanks so much! : - D I very much wanted to tell that story. It was fascinating to see different proposals -- even a sparse BLAS! -- in ACM SIGNUM in the 1970's. I'm reminded what we lost by going from an informal but archived newsletter like ACM SIGNUM, to random e-mail lists and forums that often can't be publicly searched and no longer exist.

[@literallynull]

0:00 Bro really just said to touch the grass. Bjarne is the chaddest among the programming language creators.

[@combatcorgiofficial]

The signature for this is user-unfriendly enough that I have no reason not to just stick with Eigen lol

[@abelltube]

Who is this standard supposed to help?

[@ruffianeo3418]

This is, what happens if "experts" try to define an API. It is prone to become some complicated mess, because of the experts pride. "I am an expert of REALLY complicated stuff!". Instead, EXPERT USERS should define the API and the experts should then review it, point out caveats, iterate with the expert user a few rounds and then just implement it. This has gone wrong a couple of times in the c++ standard library. Most notably the chicken scratch of . Another, non standard library examples is: Vulkan, where you need 1000 lines of boiler plate code to render your first triangle. In my decades of working "in the field", it happened on a regular basis, that people finally just threw away "expert defined" code and replaced it with something intuitively usable, written by expert users. Especially in C++ (but also in other languages), an API design pattern is needed, which does not complicate the API because of configuration options or edge cases. There has to be a way, to have simple calls, which use defaults, useful to the majority of use cases and which support "simple problems - simple solutions". The mechanism to then support the edge cases should happen before the simple calls, setting other defaults for that specific edge case. In contrast to packing all the configuration options into a cryptic maze of template arguments. This is really a comment about API design and how C++ STL usually gets it wrong. I acknowledge, that it will be harder to provide an API in a language, which lacks features, such as compiler-automated loop fusion. This should be stated by people who come up with complicated stuff, just to achieve API level loop fusion, on top of an ignorant compiler.

[@vedantranade1815]

also the fact that this is C++26 , which means most people wont use it before 2027+ , is kinda weirding me out as well. One thing which I would have liked to see, is a library which gives readable error messages, but this talk doesnt talk about it at all.

[@jonathanberkeley4109]

@@vedantranade1815 C++ is the wrong place to be if you like readable error messages

[@tiranito2834]

@@abelltube absolutely noone tbh... I mean, just look at their function names.... they avoid operator overloading to prevent ambiguity, and here they come, 0 fcks given, and introduce "std::matrix_multiply" lmao, I can't even bro.

[@bernhardmanfredgruber7156]

I think the biggest challenge for this proposal is outside the standard. It's nice to have a standardized interface to call linalg routines, but how does the implementation (i.e. the compiler + stdlib) choose a "backend" BLAS? Like, how do I tell g++ which backend to use? Will it come with a default BLAS implementation? How can I make g++ use my custom vendor BLAS which did not exist when my g++ was released? I am worried we may see a similar situation to std::regex, where the standard implementation is terribly slow compared to third-party packages. Also, if a compiler chooses to depend on an open source implementation, we may find a similar situation to g++ merging a TBB backend for the parallel STL. That was a big hassle when Intel broke some interfaces while moving to oneTBB, and libstdc++'s parallel STL, the standard library!, was broken for some time when having newer TBB versions installed. Also: how will this be exposed in CMake? I am currious (and a bit worried) what compiler and tooling vendors will do :)

[@AlfredoCorrea]

Well, I think at most there will be linker errors and underscore vs no underscore decorators missing. but other than that it is unlikely that the BLAS backend will break as TTB does. The main challenge I see is that the BLAS backend is not complete (eg, transpose complex matrix-matrix), and strides are constrained to be 1 for matrices. There will be a performance mismatch if these cases are called that STL will have to implement “manually”.

[@janciesko7994]

Great talk!

[@MikaelMurstam]

will these be faster than GLM and will they work as seamlessly with OpenGL?

[@luisclementemesquita5511]

I remember using Linear Algebra libraries in Fortran in the days of JCL, punch cards and IBM 1620. In 2024, C++ needs to hurry up and catch up.

[@franciscogerardohernandezr4788]

My wishful thinking would be to include Blaze in the standard. IMHE, it is very intuitive to use and its author Klaus Iglberger has done some serious template work to make it a great option.

[@QuintinMassey]

Are there any foreseeable drawbacks to bringing linear algebra into the STL? I don’t know enough to have my own opinion. I get that it could help with easing requirements for building your code against 3rd party linalg libs.

[@wusuoliu5431]

Baby finally!

[@sirhenrystalwart8303]

column major vs row major is just a convention, provided you store the entire 2d matrix as a single contiguous array which is the case in the first example. Nothing stops you from interpreting a C array as column major if you wish. In fact, Eigen does exactly that by default. It also seem rather disingenuous to say that this is the status quo we need to improve on when libraries like Eigen, Armadillo and Blaze are in more common usage than raw calls to BLAS.

[@FreeVFX]

It's honestly very difficult to see this being better than Eigen. Either in vector optimisation terms or library design terms. My hot take: Without an attempt at opinionated operator overloading, this just becomes an obfuscated (metaprogrammed) version of a C linalg library.

[@MikaelMurstam]

it doesn't look like you are using operator overloading?

[@Alexander_Sannikov]

"if your problem is small, then you probably don't want to be calling the blas library" well excuse me, if my problem is NOT small, I need to be using my GPU instead while designing the entire program around this fact. this leaves a conundrum: this is a proposal for functionality of the standard library, yet I fail to see any actual usecase where using it is optimal today.

[@modernsolutions6631]

There are many codes with matrices which can't just be done on GPUs. Those still need fast matrix multiplies. GPUs aren't good at every use matrices have for solving real world problems. 😂

[@MrTargetSan]

Simple question. What would happen if BLAS bumps major version due to some major improvement which requires breaking changes in API? C++ won't be able to do version bump.

[@r-jwolthuis7243]

Sure... If a different BLAS version, or any other implementation for that matter, is available to your vendor, they can make the interface matches inside their standard library. Choosing the optimal implementation for your data. And optionally give you access tot specific implementations using execution policies (like for the GPU offloading example). A different API for an external library won't effect this proposal, just it's implementation.

[@niklkelbon3662]

BLAS in standard library is good idea, but i have few questions: Main thing about BLAS is performance, but as i understand you can only pass mdspan into it (users types not supperted), its maybe ok because mdspan is just poitner + sizes But what about this 'scaled(...), transposed(...), scaled(...)' thing? I dont think its a good idea to make it 'views -like' way, if its lazy, then matrix multiplication will work with?.. What? Not contiguous memory, bad memory layout etc Its better to make O(N) swaps to transpose + O(N) scale (multiplication) than do it lazility O(N^2.71) (different N) times in multiplication? Where im wrong? Library should optimize those operations after creating, but i dont think it is possible in current interface and i dont think standard libraries will really do it, if someone will not just write it in libc++

[@bernhardmanfredgruber7156]

AFAIK the 'scaled(...), transposed(...), scaled(...)' turn into accessors (or potentially layouts) for the passed mspan. Inside a linalg implementation, you can just look at the type of your mdspan, detect e.g. a scaling accessor, retrieve the scaling factor and unwrap the mdspan. That would leave you with raw a pointer, sizes, and a scaling factor in the end. The 'scaled(...), transposed(...), scaled(...)' can thus be stripped again before dispatching to a BLAS implementation.

[@AlfredoCorrea]

the practical goal is to use as little as possible number of calls (e.g. 1) to BLAS functions with as little (or none) allocations. These “views” make sense because the corresponding BLAS functions take flags that recognize these lazy/implicit operations as part of their fundamental operations. For example, for good reason, there no function in BLAS to just multiply two arrays, the actual function is to “multiply two arrays, transposed or not, conjugated or not, scaled by something”.

[@eliasboegel]

You may find some answers in Christian Trott's talk on mdspan.

[@Spielix]

My understanding is that the implementation should be able to find out if it can call and how it should call the vendor libraries by inspecting the mdspans. Not sure how it works exactly, but only because there is e.g. a scaled accessor that doesn't mean that the scaling factor can't be extracted from it and passed to a vendor BLAS library. If you use a custom accessor instead of "scaled", it certainly can't give that to the vendor BLAS libraries.

[@cedricholocher7570]

By default it seems to be lazy. But the user can manually copy to contiguous memory. It seems the BLAS will only be called if its contiguous so that might be a useful thing to do.

[@headlibrarian1996]

linalg? Do we still need to save bytes with obscure abbreviations?

[@franciscogerardohernandezr4788]

You need a professional service to translate your humble linear algebra expression into working code.🤣

[@sirhenrystalwart8303]

linalg is hardly obscure. Much easier to read and type than the excessively verbose linear_algebra.

[@Spielix]

I would agree if it were called "la" which is what some mathematicians may name it. "linalg" seems nicely searchable for me and "std::linear_algebra" would get flamed from other people like "std::ranges::views".

[@destiny_02]

i say obscure short names are still better than verbose monstrosity like chrono

[@DKFox-rl6jq]

Which will be the first compiler to support this? Also, avoid using those underscores.

[@Alexander_Sannikov]

24:50 doing matrix multiplication on the GPU as an execution policy of a function that multiplies 2 CPU side matrices is simply not how you do it. If you want your code to run on the GPU, you can't just store its data as a span in RAM. It has to be a much deeper level of abstraction of the storage itself which is conpletely unnecessary for a standard library.

[@destiny_02]

Its not in standard, just nvidia compiler specific.

[@MarkHoemmen]

Thanks for your interest in the talk! I would encourage you to try the experimental version of this library in NVIDIA's HPC SDK. It wraps cuBLAS and therefore is able to run on the GPU, using matrices whose memory lives on the GPU.

[@imoe0x63]

when this is coming, is it in c++23 or beyond?

[@bernhardmanfredgruber7156]

Targeting C++26.

[@CharlesHogg]

C++26

[@valizeth4073]

C++23 was "finalized" months ago

[@johanekekrantz7325]

Saw this and got so excited, unfortunately my enthusiasm basically died when i saw the clumsy blas style interface instead of the more modern interfaces with operator overloading and expression templates (like Eigen and Blaze etc). BLAS is just so clumsy, verbose and far from normal mathematical notation. Its not by accident that the newer libraries are taking over in many fields.

[@vedantranade1815]

exactly! 15 minutes in I was waiting for it to end!

[@xbz24]

nice

[@robertolin4568]

Don’t you think all it takes is a good enough package management standard to include existing implementations, instead of putting all of them into the standard?

[@baka_geddy]

We need std::net for C++26

[@bocckoka]

of course it does, there is no package manager.

[@theDemong0d]

I don't really think this makes sense. Put GLM, DXMath, and Eigen side-by-side and you will immediately see how varied the requirements of linear algebra libraries can be. Trying to standardize it seems like a bad idea.

[@What-he5pr]

Why?

[@akbulutm4]

I had to scroll down too much to find this comment.

[@yellowajah]

Because they haven't gotten to the part where Jeff Goldblum sanity checks the decision making process...

[@realisticlevel2553]

I was an intern at one of the most powerful supercomputers in Europe, and when I got there and saw that they used *a lot* of fortran I was like: what? now I understand that BLAS was (is) really tightly coupled to architecture also: Why this super-generic syntax? why not rely on Eigen that is already there on C++...??

[@OhhCats]

long overdue.

[@yorailevi6747]

Unless this is going to be deprecating blas libraries this will mostly be useless

[@eliasboegel]

The intention for std::linalg is more to provide a unified interface to BLAS built for use with mdspan

[@Spielix]

How and why would the C++ committee deprecate vendor libraries? The whole point is that you can write portable code against this interface and let e.g. a vendor compiler like Nvidia's nvc++ do the work of calling the right BLAS library under the hood s.t. you only need one version to get decent performance on any platform from accelerators to laptops.

[@75yado]

Somehow I do not feel the need for this as LEDA, CGAL and other libs already exist for decades

[@johnybaker9878]

C++ OS when?

[@superuser8636]

Why can’t we just call the scalar an eigenvalue right away. Great stuff but it blows my mind this already isn’t in C++

[@Alexander_Sannikov]

at 25:00 i still don't understand why that matrix product computes alpha*A*Bt+b*C and not alpha*A*Bt*b*C

[@modernsolutions6631]

becuase dgemm is about being one optimized kernel you can use for a lot stuff your variant can be expressed as two dgemm calls while your variant can't express one dgemm call. dgemm is more versatile than your variant. Also in your variant alpha and b are doing the same thing. That's a waste

[@XEQUTE]

std::linealg 😺 !!!

[@IgorKaratayev]

With C++ ABI problems there is nothing good in putting everything to stl

[@MonsSeptime]

C++ madness at its finest! Linear algebra is now a topic for the standard library! Like 99 out of 100 C++ users are scientists who suffers without linear algebra in STL.

[@yellowajah]

why in the ever-living goodness would you want this in the standard library? Is CPP going to add an orbital mechanics library to std? an elec-eng math library? Finite Element Analysis? Does CPP ever stop to think whether or not they should?

[@niklkelbon3662]

It must be in standard library, because its fundamental things, which will not be deprecated (its not json)

[@ChrisCarlos64]

I'm not sure what the issue is exactly. If it's "Why should they add it?" Well what are you exactly against? Adding some new header with functionality or something else? This additional header is more or less additional functionality available to people who want to start using stuff in the standard library and expands the capabilities of the STL for those who just want to code. Maybe you have a library you're familiar with and love, and that's great, but often times people don't know where or what to choose when it comes to adding a new library, and then the other part is they don't know how to build said library and it adds additional complexities to just "build it once and run" for many. This may not be an issue if YOU know what you want or need, but for students, newcomers, and people just hoping to play around with the language and learn something new with minimal friction, this helps alleviate some of those pressure points for them. That's how I look at it. It hurts nothing to have or add in, and helps the language be a little more modern for new programmers who have experience elsewhere or nowhere.

[@luciusvorenus5242]

Linear algebra is used almost everywhere, especially in the applications you mentioned. And blas (which is partially implemented by std::linalg) is a long-lived standard for it, so in some ways it's like not to have std::sort.

[@LostieTrekieTechie]

Why does C++ include floating point numbers? If someone wants to do floating-point math they could implement it themselves or use an external library. Why do they have functions to deal with strings, people should just deal with byte arrays.

[@gast128]

Fallacy argument. C++ serious lacks in standard libraries (e.g. it lacks object persistence; JSON; XML; cryptography; networking; database; graph). It hampers application development compared to .NET or JAVA where this is ready available. This is a good initiative to have at least some sort of lineair algebra ready available. If one has extreme needs one can always turn to external specialized libraries.

[@alexisandersen1392]

this is too high level for the STL imo. We have Mathematics libraries that satisfy this need already, this doesn't actually affect the language, which is what I thought the STL was concerned with. Additionally this completely misses the much more general and useful concepts that arise from geometric algebra, and using multi dimensional analysis as a type system, and then matrix array etc... these all each can be represented as a sum of these multi dimensional typed objects. You could also use these sums as a model and mechanism over composition over inheritance. etc etc... there are far more useful things to add to standard library than a bog standard implementation of a linear algebra library that we have already had for decades.

[@sporefergieboy10]

Hey Siri, pull up the source code for std::lerp

[@LeeJordanMusic]

Really? Was "#include " too much typing or something?

[@maliksmith8319]

Don't really see the point. If you need blistering fast linear algebra solution you most likely know already where to get it and how to use it.

[@OhhCats]

I don't like to use external libraries. I try to use standard library when it is possible - such approach makes projects simpler. In many cases you don't need "blistering fast", you just need "fast enough".

[@benloud8740]

Whyyy. What next. A scientific library? Optimisation? Machine learning? Computational geometry? Where does it stop. Not everything needs to be standardised. Libraries are fine for stuff like this. Standards are needed to abstract away the OS for portability. For everything else, libraries are fine. Keep core C++ small please.

[@modernsolutions6631]

Becuase Julia has it

[@ferasboulala6220]

std::llm 😂

[@martingeorgiev999]

They also need to add std::calculus, I am sure it will be very useful 🤡

[@drmikethom]

What else are they going to add to C++??? It has grown to a complex dysmorphus beast. BTW, matrix multiplication is NOT pertinent as the example to motivate enlarging C++ with standard algorithms. H/W optimized Matrix Multiplication routines exist for all known platforms. The question is how to parallelize more irregular or sparser algorithms.

[@Spielix]

I think the people working on this would argue that exactly because "H/W optimized Matrix Multiplication routines exist for all known platforms" there is a interest/need for a standardized C++ way of calling all these libraries. For irregular/sparse algorithms it is harder to standardize something as there are so many different algorithms/techniques for different sparsity patterns, formats etc. and different vendors might implement different subsets.

[@MrTargetSan]

God please don't! Another library rotting in STL swamp because of standartization

[@simonl1938]

Overly complicated and user unfriendly, the usual C++ style.

[@0wO]

Whoever is in charge of designing the C++ standard is a huge failure The standard library is so big and bloated, and most of it is *incredibly slow* Keep it simple stupid, linear algebra is GLM's and other independent libraries' job

[@gast128]

A step in the right direction to have at least some sort of linear algebra readily available. Perhaps it can even define matrices / vector interface for other libraries (e.g. OpenCV) so that you don't have to learn that again for every used external library.