Matrix Computation in Wolfram Language

So welcome everyone. I'm going to be talking about matrix computation in the Wolfram language, which is a very, very big topic to try to stick into 20 minutes, which is, you know, presented a bit of a challenge on exactly what to say. I decided to start off with a bit of a historical survey. And this shouldn't be viewed as a complete or comprehensive history of matrix computation in the Wolfram language, but just sort of the parts that are either interesting to me or part of the story that I'm going to tell you today. So matrix computation and linear algebra have a very long history going back all the way back to version one in 1988 when we had things like pseudo inverse, an inverse that minors transpose, your basic solving functions, linear solve, null space, row reduce, and the Eigen family. Then in version 2 we added some more decompositions and the matrix exponential, some more decompositions in version three. Then version four was a very, very important release because we introduced PACT arrays. So, you know, one of the things that set Mathematica apart from competing systems, at least, you know, I was in elementary school when this happened, but as I understand the history, you know, Mathematica could do linear algebra on sort of anything. Symbolic matrices, mixed mathematic matrix types. It wasn't like lapack, which was sort of the, and in some ways still is the standard where you have to coerce everything to a fixed type and select the correct matrix function for your particular algorithm and data type. But this generality also made it slow. So in version 4 we introduced packed arrays, which I won't say too much about, but if you had completely uniform data, all machine sized integers, machine sized reals, machine sized complexes, then instead of using the general expression structure, we could transparently transform it to the sort of the raw machine array and then we could do linear algebra much faster on it. And so basically we have functions which if I have time I might chat about a bit to control whether you're using a packed array or not, but also transparently the linear algebra functions, they took your matrix, they would try to pack it. And if that worked, they could do the linear algebra operations much faster by working on the raw data and returning that. Another interesting thing we introduced in version 4 was to trace, which seems like a very simple function to come so late. And I will have a bit more to say about it in my next slide, I believe, but I will just cite this as the first convenience function. This wasn't a hard thing to compute on your own, but it's still nice to have a simple function for computing that. So version 5, we added the singular value decomposition, the full decomposition, as opposed to the earlier singular values, and we added more convenience functions. Another big innovation in version 5 was sparse array, which was sort of our first sort of specific matrix type with a specific data structure underneath it. Then version 6.7 and 8. I don't mean to disrespect any of them, but there were lots and lots of smaller functions that matrix predicates or generators like rotation matrix, also matrix plot, a nice way to visualize things, and Kronecker product. Then in version 9, we introduced matrix function, but more importantly for my story, we introduced the structured arrays or our first structured array, namely symmetrized array. And I'm going to end the little history lesson there because we now sort of have all the pieces we need for my story and sort of. The two themes I want to stress is that, you know, Wolfram language makes linear algebra accessible to a very large group of people, but also that exploiting type and structure can be extremely powerful and this can be done automatically or manually. So I already mentioned in the case of packed arrays, you can manually pack things, but also many functions will automatically try to pack data for you and you use that. On the other hand, if you want to use a sparse representation, you need to manually construct a sparse array and then you can do various things on that. I'll talk some more about accessibility and then I will talk about some recent features in matrix computation and then I will talk about at least one possible future direction that we're actively pursuing, so accessibility. So as I already mentioned, just being able to perform linear algebra and sort of an arbitrary matrix was a single unified interface was a great advance in 1988. Over time we added not only more powerful functions, but also simple functions that just made common things easier. I mentioned that I would have more to say about Trace. Here's a table from. I deleted a few rows to concentrate on what's needed on the core idea. But this is a table that's in our current documentation about, you know, common operations on matrices, transpose, inverse, determinant, and trace. And this, this table traces its lineage all the way to the very first edition of the Mathematica book. And this is how that table appeared in the first edition of the book. Besides the retro Courier font, which I think is fantastic, you see that the trace row was just replaced by some code for summing the diagonal elements. It's not that we didn't know that people needed to compute the trace, but the mindset in the early days was why have a function for something that's so trivial to do by hand? And in the era where you might be lucky to have four or eight megabytes of ram, not making Mathematica even larger than it already was was sort of a very valid concern. On version four, that's 1999 now we're talking about computers with maybe 64 megabytes of RAM or more, even on a high end computer. Certainly then we're like, well, let's put in the things that people expect, even if they could do it by hand really easily because it just makes life better. We've added, I mentioned a few of the others, like characteristic polynomial and others that we've added over the years. In version 13, we did a major refresh to the documentation and added a bunch of inside examples. You know, I encourage you to go to the guide page for matrices and linear Algebra, which links to all sorts of other pages that you might want to, which list more operations. But just for example, with eigenvalues, actually, eigenvalues is not a great page. Let's go to eigensystem. We have the geometry of eigensystem and so explaining what an eigenvector is, positive negative eigenvalues, or how you can understand it in terms of principal axes and the corresponding ellipsoid. Also some very nice examples to applications to quantum mechanics and coupled linear systems. Or another example of this is singular value. Whoops. Let me just do it this way. Singular. Sorry, I apparently have lost the ability to type singular value decomposition. Again. We have these examples on the geometry of, of the, of the svd, sort of what each one of the matrices does to various things or, you know, how you can understand SVD as minimizing the norms of the quadratic form. And so we've added a lot of documentation which not only tells you how to, you know, type in the command and use it, but also what does this do, why you might want to do this? And so on. Another thing we did was to create a free online course in linear algebra, which I have the link here. I won't click on it just now, but it's certainly using our tools and promoting the Wolfram language. But it's a standard sort of introductory linear algebra course at the lower division undergraduate level and sort of just putting linear algebra out there. Okay, so that's a lot about accessibility. Let's talk a little bit about functionality. So Centered Intervals is a relatively recent addition that was introduced in version 13 that represents a subset of the realtor complex numbers with the center and the bounding radius. And I put radius in quotes because it might be a complex number telling you how far you can, it's really a rectangle, how far you can go in the real direction, how much you can go in the imaginary direction. And output intervals are guaranteed to be a bound on possible answers, not necessarily tight. So if you subtract a centered interval from itself, you don't get just sort of normal zero, you get another centered interval whose center is zero. And in this case the radius is four. Because if you have a number between minus one and three and you subtract the number between minus one and three from it, then you're going to get a number somewhere between -4 and 4. And so you get that. So in this case we got a tight interval. In this case it's not so such a tight bound. And because if you look at this interval, I constructed it so it's in the upper half plane, which means we should only take the exponential of it. We should only have things in the upper half plane. But yet this bound clearly takes us into the negative half plane. Not necessarily tight, but it is guaranteed to be correct. I don't mean to cast any aspersions on this, it's just to write something that's absolutely right for every single function. I mean X maybe we could do a little bit better. But if it has to work for a bunch of different special functions, getting the exactly right bounds is a really, really hard problem. So in version 14 we're adding support to various linear algebra functions, characteristic polynomial, that eigenvector inverse linear solve and the matrix exponential. It's not a huge list, but. But it is a list. Part of why it is relatively short is that, you know, as I've already mentioned, structured is important. For example, why isn't matrix rank here? Well, sort of any invertible, any non invertible matrix can be made invertible by adding a pretty small perturbation to it. And since we have uncertainties in all of our elements, it's not necessarily a meaningful function. So there's both implementation details, but also sort of conceptual difficulties with some functions. But let's give an example of something we can do. So here's just matrix with some centered intervals and a compute can just give it to eigenfunction. Again, as you might expect, the eigenvalues have uncertainties in them. The eigenvectors have uncertainties of them. Now we should expect that we still have the similarity relationship. We get this big Sort of ugly looking matrix, although it does look like the entries are pretty small. In fact, if you ask, is zero an element of every single entry in this matrix? Well, the answer is yes. And since I'm a little low on time, I'm going to skip this little example. But I just showed how the bound on exp above was tight on some sides, but not the others. So finite fields were introduced even more recently in version 13.3. And again we're adding support for some linear algebra functions. We get characteristic, we get a bunch of very core linear algebra functions. And if we just compare the list and we've gained matrix rate, null space, row space, these are all exact things. So this makes sense. We've lost eigensystem, which I think is in principle doable, but probably hard, and matrix exponential. When you're working on finite fields, the analytic matrix exponential doesn't necessarily make the most sense. So again, let's take a look at a concrete example. Here's a matrix consisting of finite field elements. We can do the LU decomposition, which gives us a lower upper matrix, a permutation and a condition number, which is irrelevant because this is an exact matrix. We can extract the lower triangular part from the lower upper matrix. We can extract the upper triangular part and we can confirm the factorization. I'll just point out that arguably we should have actually done this and add the an identity matrix of, you know, identity elements in the finite field. I mean the identity element, the multiplicative identity, always has index one in the finite field. So one, one was as good as the other one. So those are some things we have added recently. Another thing we've added is structured matrices. So symmetrized array was introduced in version 9 and that was our first structured matrix. We then added version in version 10.1 quantity array and we over the sort of 10 cycle. We also did a lot of improvements to the framework as we worked on these two for structured arrays. And then I have to confess that the project stalled for a little while for various reasons. But we started coming back to it and we introduced our first set of structured matrices and 13.3. Each release since them since then has included more structured matrices. I'm not going to list them all here. I have a link here to the complete list. But what's new in 14 is Hermitian matrix, symmetric matrix, unitary matrix and orthogonal matrix. These are designed to be user extensible. In fact, many of these newer types, in fact, I believe all four of these types are added in a way that a user could have done it. We've done it as part of the system, but it was done in exactly the same way that a user could have done it. We don't have quite all the kinks worked out and we certainly don't have the documentation written. But that's. I believe Devendra mentioned in the earlier talk, that's something we're working on. So why do we care about structure? Well, I mentioned earlier, performance can be pretty dramatic. And to just give you one example. So diagonal matrix gives you a normal matrix for small dimensions, but a structured matrix for large dimensions. Now, you can control this using the target structure option. I chose not to display all 25 million elements of the matrix, but if I look at the head of this matrix, it is a list, it is a normal matrix. Inverting this matrix takes a little bit of time, about four seconds. But as you might recall, inverting a diagonal matrix is really easy. This should really be faster and yeah, 0.17 seconds. So, yeah, a little bit of a speed up there. Only two orders of magnitude there. And every structured matrix knows what algorithms it's been optimized for. And you can ask it. Well, and so here are the algorithms that diagonal matrix has been. We've added dedicated code. So structured matrices sort of occupy this middle road between the two extremes of sparse array and packed array in that they're sort of whether you get one can depend on your inputs, but you can control it. So it's sort of in between whether it's manual or automatic. And so those are some things that we have ready to show you guys. Then something that we're working on is GPU computing. And Tom Wickham Jones gave a talk about this earlier today, And in particular on Coda, excuse me, on CUDA. And we've long had CUDA link and OpenCL link, so you could do GPU computing. But this worked well but required a lot from users. You had to know how to manage memory yourself. We're now looking to package this to be more accessible. And the first front of that is being released in 14. And Tom talked about this CUDA vector, CUDA matrix, which sort of makes things a lot easier. Now I'm working on sort of the linear algebra side. And so what I'm going to be showing you is results of some performance comparisons of hooking up linear solve to the gpu. And I want to be very clear that this is an early experiment. This version of linear solve ran on precisely one computer in the world because I literally hard coded the. The paths to various libraries as you know, as they were on the file system of that one computer. Because like I said, this was an early experiment. I'm working very hard on this. I hope to have something to show soon, but maybe for the spring, but I can't promise that. And so we took a packed array so a bunch of machine reels and we perform two operations at double precision on both the CPU and gpu. So I was using a very high end GPU in RTX 4090 which costs I think like 16 or $1700 and high end CPU, a dual cores, latest generation dual core Xeon. And you can see that for lar and I basically just did linear solve on matrices and sizes of 1024 up to 46 times 1024. You can see that for large matrix sizes we got some speed up, but it was nothing dramatic. In fact, on lower end cards the GPU never beat the CPU because that's really not what they're great at. Then we did a different test where we took the same data. We did the CPU computation as normal for the GPU computation and this time is included in what's being displayed here. I took the time to convert that to single precision prior to sending to the gpu, then converting back to double precision to return the same data type. Obviously you're going to lose some accuracy, but for many applications you don't need all those extra digits. So if you don't, well, you can get some nice split up. And I don't know exactly what's happening here, although I suspect it has to do with whether everything fit in the card memory at once or it needed to batch the computation. But you can still see you get a very healthy performance speed up. Even on lower end cards it's not quite as dramatic, but you still get a very noticeable speed up and that's what these cards are very good at. So, you know, future direction, something we're working very hard on. There are all sorts of logistical issues to work through, but that's what we're working on. So in conclusion, I just want to say that, you know, Wolfram language has a long history of bringing innovation to matrix computation and linear algebra. And that's both greater performance and greater accessibility and we're continuing to innovate with that. I thank you for your attention and I'm going to try to figure out how to look at the Q and A. Okay, So Anton has. Comment about he has a package for sparse matrices with named rows and columns. I believe that's certainly something that we're working on and maybe Rob Knapp has is giving a talk on that. I'm not entirely sure if anybody knows. Feel free to comment in the, in the question, in the Q and A section, will there be generalized functionality to compute eigensystems for tensors? So yeah, that's the. So unlike the rank 2 arrays, where there is a single well defined notion that everybody agrees upon for eigensystems of matrices for higher rank arrays, there are multiple different notions and they don't have quite the same. I mean it's very important, but they're different notions and they all have slightly different meanings. It's not quite as straightforward as it is for matrices. You know, that's certainly something we're interested in, but it's not something we're going to be delivering in the very short term, as far as I'm aware of. So there's a question about packed arrays in version 4, when the compiler project started only in the last few years. So yeah, so let me just briefly mention. So if you. I don't know, if you just search for pactray, if this will come up. It will come up. We had functions in the developer package to convert things to and from. Boy, there are no examples on this page, which is a little disappointing. Because I need. All right, let me. If I do two packed array of 1, 2, 3. Well, that looks just like a normal matrix. But if I asked is this a packed array? The answer is yes. Packed arrays looked just like normal matrices and unless you use these developer package functions, you couldn't necessarily tell the difference between them. There's also, you could sort of, you could use, you could do messages to see when packed arrays were automatically unpacked. And if you had a very larger array, you know, if you unpacked it, it's going to be using much more memory, at least the factor of eight. So it was very useful to be aware of where those things with the compiler project, we're sort of exposing that as a definite distinct type that you can compile for. And it's basically the same packed array, only it's now sort of exposed more directly in the language, whereas previously it was hidden and let's see, With named rows and columns. Yes, that's certainly something we're working on. Right. So I, you know, I don't have a. There's a question about intuition regarding matrix computation with finite fields. I don't have something prepared beyond what I already wrote. So I will refer you to the earlier talk by Adam about finite fields in general. But when you're using finite fields for coding, as in, you know, compressing data or cryptography, that kind of coding. That's when you might be doing linear algebra with finite fields. And then I guess there are requests for identity matrix of zero, you know, Einstein summation. Yeah, those are. Those are. Yeah. Well, I would dispute whether all of those are matrix computation and not general tensorial operations. And there are all sorts of dragons. Those are not simple. I'm happy to have an offline conversation about those, but those will be the aptly named topics not appearing in this talk. To quote the Holy Grail. So I think those are all of the questions. I will give people a another few seconds to ask a question. Is there a predicate function to test if a system of equations is linear? I don't believe we have a documented function that does that. We have some internal functions that do undocumented internal functions. If. But I don't believe. I don't believe we have like linear queue or. Yeah, we. Yeah. So if we had this function, it would be here and it's not. We don't have a dedicated function for that. But that's maybe something we could expose. I'm a couple of minutes over time already and I don't see any additional questions, so thank you very much and have a nice rest of your conference.