So, I just found out three seconds ago that I'm introducing my colleague John Urschel, but John's very easy to introduce. He's a man after my own heart. He likes to do mathematics and he likes to do it with Julia. So, you know, what's not to love? So, John is a graduate from a school in the state. He's a graduate from Penn State as an undergraduate and was also a graduate student and earned his PhD at MIT. And I feel lucky that he was in some of my classes back then and I feel lucky that he's my colleague right now. So some people think that one of the interesting parts of John's biography is that he was a professional football player for the Baltimore Ravens. I don't know how to play football. Nobody else at MIT knows how to play football. They don't even know what the game is. So nobody seems to be aware of this or even care. But other people seem to think this is interesting. So, without further ado, let me introduce you to John Urschel. There he is.
All right, great. Okay, that was a very nice introduction. I'm not sure I've ever had one quite like it. Okay, again, and I never will again. So, what I'd like to talk to you all about today is an interesting math problem, a pretty cool math problem that has a pretty long history going back all the way to sort of, you know, the start of modern computers. And I'd like to talk about the history of the problem, some recent progress, and actually how, you know, Julia has really helped with a lot of that recent progress. Okay, so, Alan did a great job introducing me, but I'll say a couple things about myself. I'm an assistant professor at MIT math and really as far as research goes, I just really like matrices and anytime there's a matrix somewhere, you'll probably find me there. So numerical linear algebra, matrix analysis, spectral graph theory and I'm also a Julia user since 2016 when I took random matrix theory with Alan and my sort of approach to Julia may be sort of quite different than many of you. I've enjoyed going to lots of talks and learning about different things. I'm sort of a greedy Julia user. I care about doing math. I want to solve problems and I find Julia is an amazing tool to help me do that. It's a language where, you know, coding really feels like, you know, close to mathematics and there's something I want to be able to do, I just sort of Google around and what do you know? There's a package for that. There's a package for some random matrix theory thing I want. I wanted to do something over GF2. There's a package for that. And so, if you've ever, you know, helped create or maintain one of these packages, I thank you. Now, to the problem I want to talk to you all about today.
Well, it's a math problem, and I thought to warm us up, why don't we warm up with a math problem kind of like high school, college style. So, this might feel a little bit like the type of math problem you might have seen in high school or college. So, here we go. So imagine you're making rice and so there's, you know, there's good grain, there's okay grain and not so good grain and you're taking like combinations of the three of them. And so the problem tells us if we take three bundles of the really good stuff, two bundles of the normal stuff and one of the bad stuff, we're getting 39 — I don't know the units, but we're getting 39 something of grain. And then we do this two different ways and we get different amounts of grain. And really what we'd like to know is if I just use the really good stuff, how much grain am I going to get out? Or if I just use the normal stuff, how much grain do I get out? If I just use the bad stuff, how much grain am I getting out? And to, you know, some of you hopefully, you know, if you sort of make some assumptions about how, you know, the grain is made, you might think of this kind of like a linear algebra problem. And you could represent this using the matrix on the right hand side where the first three columns represent our three variables. The last column is representing sort of the amount of grain we're making. And maybe we want to solve these three equations for these three variables to see how much each one would make. And uh...
Yes, yes, the left side is very old Python. No. So the left side is actually representing what's on the right with sticks. So one vertical stick is one, two vertical sticks is two. A sideways stick is a five. And if you have two to the left and to the right, it's different numbers. So this leftmost column reads 1, 2, 3, 26. So that's my bottom row. Middle column is 2, 3, 1, 3, 4. And then the last column is 3, 2, 1, 39. And so I'm representing these numbers with sticks, which might hint at how old this problem might be. It's very old Python. Okay. So how might you solve this? Well, in school, I think many of us learned this algorithm called the elimination algorithm. We take multiples of one row and subtract it from the other rows to remove values and create zeros. And we do this over and over again. So for instance here what you'll see is happening is the second row is being scaled by three and then we're subtracting two times the first row from the second to create a zero in that leftmost middle entry. And again this is being done on the left with sticks.
And you can do this again to now put a zero in the bottom left corner. And then you can do this again to put a zero in the sort of bottom left middle corner there. And now with something like this you can go back and now you can solve for the variables successively. And the answer you get out, in case I know you all were dying to know, is that one bundle of the really good stuff yields nine and a quarter of rice. The normal stuff yields four and a quarter and the not so good stuff is only 2 and 3/4. Okay. So this is an example of Gaussian elimination in a very early form. So this example actually appeared in the Nine Chapters on the Mathematical Art dated about 200 BCE. So like I said, very, very old Python. And really this, it has some differences. The example I showed you is, it's pretty suggestive because no negative numbers show up. Nothing strange happens. The method can look a little different than the modern elimination method we learn in school but really in spirit it's about the same and yeah, so this was many, many years before Gauss was even born. Okay, so now let's talk a little bit about what do I mean by the elimination algorithm or Gaussian elimination or an LU factorization. I want us to think back in our brains what I mean is I want to factor some n by n matrix A as the product of two matrices L and U, where L is lower triangular and U is upper triangular. And the way I can do this is by at each step taking my first row and subtracting multiples of it from all the other rows to create zeros in the top, in the sort of first column. And this can be represented through a rank one update which you'll sort of — can you see my cursor? No. Which you'll see as the b minus x y transpose over alpha. And the idea is we do this over and over again. Eventually we get a lower triangular times an upper triangular. We solve this system and this allows us to solve any arbitrary linear system. Okay. And not every matrix has an LU factorization. For instance, if you have a zero in the top left, that could be a problem. You might need to permute your rows. You might need to permute your columns. And the issues that show up aren't just a zero in your top left entry. But really, when you have some entry in your top left that's pretty small. So, you could imagine, you know, having a zero there is a problem for a computer because you're trying to divide by zero. But in fact dividing by something small can be just as dangerous. So imagine, you know, you have an epsilon in your top left, you're going to get some LU factorization where the entries start looking like one over epsilon. And now I know I'm in front of an audience where I don't need to convince you that this is scary, but okay, I'm a mathematician. I'm in a department with lots of pure mathematicians. And the moment I say, you know, I'm thinking about error in Gaussian elimination, I get stopped immediately. They say, 'Hold on, hold on. Is Gaussian elimination not an exact algorithm? You're just adding, subtracting, multiplying, dividing. Gaussian elimination is exact in exact arithmetic.' So I know you all know better, but please just like act like you don't. So why don't we compare? So let's open up Julia. We'll use linear algebra. And I'll take a 1000 by 1000 matrix and just a length 1000 vector. Okay. And let's see how exact arithmetic versus floating point arithmetic does. Okay. So first I'm just going to run it once. I'm not going to tell you how long it took. I want to avoid precompile time. And now how long do you think it's going to take in Julia to solve a 1000 by 1000 linear system? Any guesses? I did it on this laptop. Let's see. It took — yeah, pretty close. Not too far. So it's very, very fast. That's the easy thing to know. Now using rational BigInt. How long is it going to take me? It would take me about the time of this talk. It would take me about 40 minutes. And well, how big is one entry of the exact answer? There you go. The division symbol is right around there, if you're looking for it. And then you might say, well, okay, it took me 40 minutes, but I got the exact answer. You know, how much better should I feel about myself for having the exact answer? And sadly, the answer is not that much better. Like the thing you compute is accurate to, you know, 10 to the minus 3. So somehow exactly computing things for Gaussian elimination is not really a strategy for success. So we really need to do things in floating point arithmetic. And in floating point arithmetic we really understand how error propagates in Gaussian elimination pretty darn well. And this has been known since, you know, since Wilkinson in the 60s, you know, one of the sort of godfathers of numerical linear algebra. And what he figured out — okay, this notation looks a touch complicated, but the idea is I want to solve Ax = b using an LU factorization. The answer that I get is a solution to almost the right problem. If this big matrix E I'm writing down was zero, I'd have the exact answer. How far off am I? How big is this E? Well, it scales like n, the dimension of my problem, times u, the unit roundoff in double precision, about 10 to the minus 16, and then it scales by the size of the entries in my L and U. And so really what this means is that for an LU factorization, if I compute L times U and those numbers don't get too big, I'm happy and I know that I've got a pretty good answer I can trust. So really the game when it comes to this elimination algorithm is to just say how can I rearrange the rows and rearrange the columns of my matrix so that when I do this LU factorization, entry stays small, entry stays small, I'm happy, I know error is small, reasonable enough. Okay. So how do people actually choose to permute their rows and permute their columns? So there's three popular strategies. I'd say the most popular. The first is complete pivoting. This is to say the pivot that I choose, I want it to be the biggest magnitude entry in my entire matrix. Rook pivoting sort of says no, I just want it to be the biggest entry in the first column and the first row. And partial pivoting says I just want it to be the biggest thing in the first column. And so for context, if you look at LAPACK, the standard is partial pivoting. So this is what the backslash command calls in Julia, in MATLAB. This is what linalg.solve calls in Python. So you're using Gaussian elimination with partial pivoting. Gaussian elimination with complete pivoting is sort of a famous, more stable, sort of more theoretically popular algorithm. This was actually von Neumann's like favorite. This was actually the standard technique around the 40s. And then rook pivoting is pretty modern. The idea is it's supposed to be some Goldilocks between these two strategies. Okay. And for this talk, I'm going to focus a lot on partial pivoting and complete pivoting and sort of the parallel stories that have shown up in these two areas over the years. Okay. So now we know what Gaussian elimination is. We know...
About some of the choices and strategies people use. Okay, I showed you an example in Julia, but computers haven't really shown up so much. I said it's 80 years of computing, and okay, I'm not counting 200 BC as, you know, computing, as close as that looks to Python. What I'm going to consider computing is what happened around the late 1940s.
And so at the Institute for Advanced Study, John von Neumann and Herman Goldstein, they created what was called the IAS machine. And this was one of the first stored-program digital computers.
Which one was the first? Which one was the second? This seems like a very sensitive topic. IAS claims it was actually the first. Other people claim this is not the case, but it was one of the first and very early on. And so they have this machine, and what they want to do is understand: can we do anything with this thing? Can we compute with it? Can we give it a problem and trust the answer we get out?
So they really wanted to study rounding and computation. And the very first algorithm they wanted to understand was, gosh, elimination with complete pivoting.
Unfortunately, they did not program in Julia in 1947, but hopefully, you know, some of their acolytes have learned their lesson. And okay, von Neumann, for those of you who sort of know about, you know, this famous mathematician called the Martian, unfortunately he died quite young, but Goldstein, in his sort of part memoir, part, you know, book about computing, said that they chose this topic intentionally as the very first paper they wrote to really understand how to do computations on a machine because they felt that solving a linear system was absolutely fundamental.
And in the paper they actually failed in their goal. They managed to show that Gaussian elimination for a symmetric positive definite matrix under complete pivoting was numerically stable. So you can think like some version of Cholesky is numerically stable. And in the paper they felt that Gaussian elimination in general should be stable, but they felt that proving this was beyond the tools that they had at the time. This is roughly what they said in the paper. Okay.
So then around the 1960s, Jim Wilkinson comes along, and okay, Wilkinson, a Turing Award winner, really revolutionized numerical linear algebra and, you know, parts of numerical analysis as well. And in the 1960s he produced error estimates for Gaussian elimination that look like the one I showed you before. He was the first person to recognize that this thing is actually stable as long as entries don't get too big, and he recognized the dependence on what we'll call the growth factor, which I wrote before, but you should think about as the largest entry in L or U.
For partial pivoting, the algorithm which, you know, is implemented and we all use, he showed that this quantity can be exponentially bad in the worst case, and that this bound is tight. For complete pivoting, he produced an upper bound using a pretty clever technique that says we're doing much better. It's n to the log n, which still gets big, but gets big gracefully. But he didn't believe this. He thought, no, this is absolute nonsense. This algorithm should be completely stable. The entries of L and U should never get bigger. Well, for complete pivoting, the entries of L never get bigger than one. The entries of U should never get bigger than, let's say, n for an n by n matrix. Okay.
And to illustrate what can go wrong in partial pivoting for this sort of catastrophic matrix that Wilkinson was describing. Okay, let me use linear algebra again. I'll just make n 100. And now what I'll do is I'll just define the matrix that Wilkinson was talking about. So to show you what it looks like, it has ones on the diagonal, minus ones below the diagonal, and ones in the rightmost column. And when you do Gaussian elimination with partial pivoting to this matrix, what happens is the entries in the rightmost column double at every step while the rest of the matrix stays the same.
And so that means the entries in the rightmost column are getting exponentially big, and bad things start happening. So, for instance, if I just take some x, some, you know, random vector, I want to choose x so I know what the actual answer is. I compute my b, and now I just want to get my x back and see how well I do with the backslash command.
I'm getting, you know, this is not a very good answer. Just outputting zero is almost equally a good guess. So somehow something has gone really wrong here, even for a 100 by 100 matrix. For people who, you know, are really sort of like numerical, you know, really initiated into numerical analysis, this is around the time I get a question around, well, okay, but how ill-conditioned is the problem? And the condition number is not bad, it's like 45, which is pretty acceptable for a 100 by 100 matrix. Okay.
So now what I'd like to talk about is really, I've set the stage about the history of this problem, what does the modern story about Gaussian elimination look like? Okay. And really the modern story takes two, you know, very different looking tracks. On one hand, you have partial pivoting where we know what's the worst thing that can happen, and the question really becomes how rare is this bad thing? There's a bad thing. It can really bite you. Is it very, very rare? Is it just a little rare? And the study and sort of research that has gone on is around understanding whether we should just not worry about this because it never happens, or to understand why this problem is actually not much of a problem at all and easily dealt with.
So average case analysis is the idea that, okay, if I take some average-looking matrix, this will never happen. The growth factor, how big the entries get, won't get too big. And then there's this idea of smooth analysis that maybe bad things can happen, but even if a bad thing happens, it's very, very brittle. Meaning that even just moving a tiny bit away from a very bad thing, all of a sudden leads to something where there's no problems.
On the other hand, for complete pivoting, we don't even know if something bad can happen. I mean, okay, certainly we know things can't get too bad based off Wilkinson's proof, but really there the search is to understand do the lower bounds and upper bounds, do these match somewhere and where do they match? The upper bound says something looks like n log n, which does get pretty big as n gets big. The lower bound says n, and these two things are embarrassingly far away. And so understanding where does the truth lie, and, you know, where in this band and how do we get there? Okay.
So I'll start with partial pivoting. So for partial pivoting, this idea of an average case analysis was really championed by Nick Trefethen in this paper with Trefethen and Schreiber, also in his book, Trefethen and Bau. They sort of numerically studied the growth of random matrices, and an observation really started to sort of occur. You know, Alan observed this as well, that for partial pivoting it seems like the largest entry of your matrix really scales like root n with high probability. I won't say the precise, you know, version of what people conjecture, but this is roughly the spirit.
And recently there's been good progress to show that this quantity is at most polynomial, but okay, the techniques of this 2024 paper, it's a very long paper, it's very technical and...
Oh no, that's not good. Ah, okay. It's keeping us on our toes. Comes right back on. Okay. Perfect.
Perfect. And here's an example of a picture that shows you that, you know, really this is the truth, and you can even see this behavior, you know, when you just run a couple lines of Julia code. This is just a histogram. I take 100,000 samples of the growth factor of rand n by n for different powers of two. And you can see, you know, when n is 256 all the way up to 2^18, really this scales like root n, like this really does seem like the truth, although this is incredibly hard to prove, I believe. Okay.
Now a different approach says, well, you know, instead of saying, you know, I take some random matrix, why don't I, why don't I try to find some Goldilocks approach? You know, thinking about worst case is far too pessimistic to say, you know, what's the most dangerous thing that can happen, but saying what happens on average also is pretty soft. Like, you know, what really does average mean? You know, okay, this is a quote of Alan's that, you know, it's a mistake to say that a random matrix looks like any old matrix because, you know, the distribution of a random matrix, I don't know what the distribution of matrices that we deal with or that Julia gets on a given day, but it certainly doesn't look like a random matrix distribution.
And so perhaps an even better way to try to understand what's going on is through an idea of smooth analysis. And smooth analysis in a nutshell says, "I'll let you pick any matrix you want. It can be as bad as you want. It can be, you know, the worst thing in the world. And what I get to do is I get to take whatever you give me. I get to just sprinkle on a tiny bit of noise." And the idea is that if I just add a little bit of noise, everything is going to be okay. And in fact, people have been able to prove that without pivoting, this...
certainly holds. And with partial pivoting, this should be even more true. Though somehow the actual pivoting strategy, the fact that you have a greedy algorithm, makes it harder to actually say anything concrete. But you know, as a numerical experiment, here's Wilkinson's matrix. I add a little bit of perturbation, this horrible thing that just hurt us. If I just add even 10 to the minus 16 level sort of perturbation to this matrix, already the largest entry drops from 10 to the 30 to 10 to the 18 for that 100 by 100 example, and it just gets better by then, and this is without pivoting. Pivoting should only help you.
But still, despite sort of being able to see exactly what should happen on a computer, despite Julia telling us this is the behavior, actually proving any of the things I'm showing you has proven incredibly difficult and incredibly elusive. So somehow we can really see the truth of what things look like, what should happen, but actually being able to say anything with math is really, really hard.
On the other hand, if I move to complete pivoting, it's the complete — well, not the complete opposite — but computers are not really telling us what's happening. So in the 60s and 70s, 80s and 90s, people were trying to prove this conjecture for small values. They were unable to actually improve the upper bound of Wilkinson in any sort of meaningful way, and so they focused on small cases. And then Nick Gould in 1991 managed to find a 13 by 13 matrix with growth factor 13.02 using LANCELOT, the optimization software.
Alan a year later found a nearby matrix with the same growth factor in exact arithmetic. But since the 90s, really not much had moved for quite a long time. The idea is we know that this conjecture that Wilkinson made in the 60s is barely not true, and we have pretty much no improvement to the upper bound.
And recently Alan and I, and also Alan and Ankit Basain and I, have gotten quite interested in this problem, and you know we tried to sort of bring math and also Julia to bear on this problem in the hope that something fruitful would happen.
And in fact, the first thing we did is we managed to show that this conjecture that Wilkinson made in the 60s, or Wilkinson and Crier, that this is actually not true for just most values of n, and they're off by at least a multiplicative constant. And in fact, the way we managed to do this was through a combination of mathematical theorems and using nonlinear optimization through JuMP and IPOPT.
And so the idea is we used JuMP calling IPOPT to try to find locally extremal matrices, and then we proved math theorems that allowed us to take these extremal matrices in floating-point arithmetic and prove that they must be genuine examples, and then to interpolate those to construct matrices that have large growth for all values of n. So this lower bound, a crucial ingredient, was really using IPOPT, which was made, you know, quite easy, at least for me, through JuMP. So I'm very appreciative for the existence of JuMP.
Then about a year later we set our sights on the upper bound. And for the upper bound I'd like to say a couple words about the idea of how we improved Wilkinson's idea. And okay, I've been to a lot of optimization talks here, so I thought why don't I show you a little bit of a flavor of what the upper bound looks like. It's an LP-type argument, and then I can tell you a little bit about the spirit of how you improve this, and then also how using JuMP and using goi and other things sort of helped come up with the right clever mathematical proof. Okay.
So for Wilkinson's proof, what he does is he bounds the products of the pivots of Gaussian elimination using Hadamard's inequality. So I'll let the pivots P1 to PN denote the entries that I'm using to remove rows at every step. And the point is the maximum growth factor, you can really think about the biggest entry in the bottom right corner of your upper triangular matrix. And so we're trying to maximize a ratio of P1, which is the pivot of the one by one matrix, divided by PN, the pivot of the original N by N matrix. And this extra constraint I have is something called Hadamard's inequality. It bounds how big the determinant of a matrix can be. If you haven't seen it, don't worry about it too much.
What you can do is you can convert this to a linear program by doing a log transformation. And amazingly, the simple structure of this LP makes it really easy to solve. So I can write this LP as a matrix inequality.
So I'll write it like this. And the amazing thing is that this matrix on the right — on the left, rather — is lower triangular. So it's easy for me to take the inverse of it. And it just so happens that when I take the inverse of it and I apply it to the vector corresponding to my objective, I actually get something non-negative. And so the point is, if I take this inverse, apply it to my objective, and I actually get something non-negative, I know it has to be the solution, and I'm done. I've solved my LP. No software is needed. You can just do this with pencil and paper. So this is very, very nice.
Now, our idea to try to do better than this is a pretty simple one. The idea is when you're doing Gaussian elimination, what are you doing? You have a matrix, and every time you're adding a rank one matrix to it, and you're doing this over and over again. And so in spirit, when you do Gaussian elimination, really only one of two things can happen. Either you get some matrix that has really big entries from adding a bunch of rank one things, or you have a matrix that has really large rank. And when I mean really large rank, I mean sort of like numerical rank. You know, the main singular values. It has many singular values that aren't too small.
One thing to notice is that you can't actually have both, because if I started with some matrix, maybe all its entries were at most one, I add a bunch of rank one matrices to it. Maybe I add, you know, log n of them or root n of them. Then when I'm done, either the entries aren't too big and it has pretty high rank. But if the entries are really big, the rank needs to be pretty small, because it's essentially only the perturbations that I added that make the entries big. And so the rank really has to be mainly those things. And so this
tradeoff allowed us to prove a sort of tighter version of Hadamard's inequality to try to come up with a better linear program. So we get some more complex linear program that looks complicated. We have a constraint for every choice of K and L, where you can imagine K is the size of your matrix and L is how many steps of Gaussian elimination you look, and you want to apply this inequality to that jump.
So we get this improved linear program from this idea, and then the big question that I found myself asking myself for quite a while is, okay, what's the optimum? Is this actually going to give us anything better than Wilkinson's original bound? You can do lots of things, and somehow mathematically you don't get an improvement out. And so the big question is, do we get an improvement, and if so, how do we prove it?
And the answer is, well, I use JuMP. I use a linear programming solver. I use one that works very well and scales up to pretty large values. And then I just look and see what constraints are tight. And so I made plots like that where the black dots are telling me the constraints that are tight, and the rest are telling me these are ones that are never really achieved. And then I just say, well, okay, well, let me just try to use those in some clever way.
And what do you know? The observation that really I should only care about the constraints on like this right part of the triangle and this line going out was exactly the inspiration that we needed to find a really slick and really clever proof. The proof uses an idea involving taking this finite-dimensional linear programming problem, converting it to a continuous one, and applying LP duality in a way that really doesn't look like LP duality. And it all looks very clever. And I have to say it feels like one of the clever proofs, you know, that I've come up with. But really the cleverness comes from, you know, from JuMP telling me this is the answer, this is the technique, this is the idea.
And from this, this sort of thing pops out. And so from here you get this improved bound. And this is the first improvement to Wilkinson's upper bound since 1960.
And okay, I could show you a little bit about how the proof goes, but I think I might spare you the details unless, you know, if you have a question you're dying to know, you can know. But the main thing I'd like to say is that, okay, where we're at now is we have this fascinating, you know, math problem that, you know, really started in the late 1940s and still we do not fully understand in either of these tracks.
And I think, you know, computing power — I think sort of like seeing what the landscape of this type of problem looks like is really a tool that can bear fruit. And so I hope you enjoyed the lesson on an old problem and I hope you enjoyed learning a little bit about Gaussian elimination. Again, I'll stop there. I'm not going to show you this part and I'll just say thank you. [Applause]
Yeah, so I stopped a touch early. Unfortunately I was invited to come to Julia. I was a little bit of a late addition. So I was only at it about two weeks ago I think and my son has a dentist appointment that I promised I'd take him to in the morning. So I am running to the airport at 5 but I'm happy to take questions till then.
So do we have a few questions? Somebody over there.
