Dynamical low-rank approximation for radiation transport
Frank, Martin
The dynamical low-rank approximation is a low-rank factorization updating technique. It leads to differential equations for factors ina decomposition of the solution, which need to be solved numerically. The dynamical low-rank method seems particularly suitable for solving kinetic equations, because in many relevant cases the effective dynamics takes place on a lower-dimensional manifold and thus the solution has low rank. In thisway, the 5-dimensional (3 space, 2 angle) radiation transport problem is reduced, both in computational cost as well as in memory footprint. We show several numerical examples.
 So thanks a lot for the invitation. I would have very much liked to be at CIRM, you know, like everybody else here, but we'll have to do and there is hope that we will be able to meet in person again in the future. And I guess all of us have realized what we're missing. So let me talk about a new thing that I've been working on or my group has been working on for a year with different collaborators who are listed here. And it's I think potentially also bigger idea for numerical methods for kinetic equations especially. So it's called dynamical low rank approximation. I will explain it using radiation transport because that's what I've been doing, but it's a general method that I believe is of interest. And I should say I'm not the only one or we're not the only ones who have realized that. So let me explain what it is because it's nice mathematics. So for the moment, imagine a very simple kinetic equation. And for the sake of simplicity, I'm going to assume that it's just 1d plus 1d. So it's 1d in space. And it's 1d in velocity or in radiation transport. It's an angle variable. So you have a probability density depending on time, one spatial variable and one angle variable. X might be an R and the angle in radiation transport is usually in the interval minus one one, but the details don't matter. And here I've written down just just just the just the 1d kinetic equation. It's 1d linear transport. So you have advection here this term. You have a scattering term here. And the fly is just the integral over velocity space averaged. There is a total cross section here loss of particles absorption scattering and there is a source here. But again, the form of the equation doesn't matter for now. I've put everything on the right hand side because I'm going to consider this to be an ordinary differential equation in a second in some functional space. Alright, now assume that that you discretize this equation somehow or semi discretize it somehow in in space and in velocity in these two one dimensional variables. And what I end up with is is a quantity that depends on two indexes, right? The index I is the spatial discretization, the index J is the angle of velocity discretization. Alright, so I I end up with a matrix. And if I don't discretize in time, then this is a matrix differential equation that I get. And in general, it can be any right hand side here for the next couple of slides. And it can be nonlinear, if you wish. Yeah, so why is your discretized kinetic density? And this is the right hand side you get. And obviously, this is this is something quite natural that you get you had these two indexes in kinetic theory. And you can also imagine that these matrices become quite large, right? So if you do a 3d plus 3d kinetic equation, you have more indexes than this, but you can put all of this into into a matrix, right? And your one direction is space and the other one is velocity. Alright, so so it's called dynamical low rank. So what are low rank matrices that's your linear algebra one or two lecture? Well, rank is most easily defined by a singular value decomposition. So that's what I've written down here. So I'm speaking of a rank r approximation. So this y gets a little index r here. And this looks like a singular value decomposition. So there is a matrix u times a matrix s times a matrix v. And an SVD usually looks like this, right? So you have your object, which should have been called y here. You have you have a matrix which looks like a towel, right? So this is m by m. This is an m by r matrix, this s is r by r. And this is r by n. If the original matrix was m by n and r is the rank. And in a singular value decomposition, this would be a diagonal matrix. But for now, we allow these matrices also to be full. So that's if I can do this now, then that defines a rank r matrix or a rank r approximation. If I take the full singular value decomposition and I just cancel out some of the singular values, then I have a rank r approximation of that matrix. Now, something very well known from your standard linear algebra course. I should mention a related concept here, which is sparsity. So rank and low rank is related to sparsity, but sparsity is a little bit different. Because in sparsity, you would ask for a matrix, for example, this s matrix here to have a lot of zeros, right, which is a different concept than asking for a low rank. So I will maybe discuss this in the end and maybe there are questions about this. Because ultimately, the goal is to save some computation, to save some memory. All right, so that's the definition of a rank r matrix. And the idea is that we want to approximate this matrix from our discretization just by a lower rank matrix, right, by a rank r approximation. And so that's what we do, right? So I view this as an ansatz, right? I basically say I assume my solution is of low rank, which means it has this SVD-like decomposition, and I insert it into the equation, into my abstract ordinary differential equation for the matrix that I have. Obviously, this truncated or approximated low rank approximation will not satisfy the equation anymore, right? If I insert it, I will end up with something that's residual, something nonzero. And now comes basically the trick of this dynamical low rank method. I do something that I also do in a Galerkin method. And I would explain it in such a way that this is now a nonlinear Galerkin projection. So in Galerkin, you would say the residual I demand is orthogonal to the ansatz space. That's the standard Galerkin thing that you do. Here, we do something different. What we do is we take the right hand side and we project it onto the tangent space of rank r matrices. Why do we do that? Well, this matrix doesn't satisfy the equation. We need to make sure that we get an equation that it satisfies, and we need to make sure that we do not run out of this space of low rank matrices, which is now not a linear space anymore, but it's a nonlinear manifold. And so in this way, I make sure that if I start with a low rank matrix, if I project the right hand side onto the tangent space, then it basically says the time derivative of this thing is in the tangent space, then I remain for all times in the space of low rank matrices. So that's why I need to project the equation. And again, this is a nonlinear object. If I'm at one point in this manifold of low rank matrices, which is y r, then at this point, there is a tangent space. And that's what I do for now. Now, the question is, what is the tangent space? Can I characterize it? And can I express this projection in such a way? Because ultimately, I want to design a numerical scheme. And the answer is yes, I can do that. And I want to give you a short argument of how the tangent space looks like. It's not a proof, but it's the essence of the proof. So if I have an SVD-like decomposition of a matrix, y is USV transpose, then you can do a simple product rule, if you wish. So if this is a variation, delta y, think of it as a derivative, then I have a product rule. I have the derivative of the first one or the variation of the first one, the second one, and the third one here. So this would be an element in my tangent space here. And then I can say something about these variations, because I want u and v to be orthogonal. So that means u transpose u is the identity. If I differentiate this, I get this equation here, namely that essentially my variation delta u times the u matrix, the original u matrix has to be zero in a symmetric way here. So this equation. And the same holds true for v, because it's also an orthogonal matrix. And for s, I just need that the delta s is also an R by R matrix. All right. So then what I do, I make this a little bit simpler. I satisfy these equations by demanding that each of these two terms in the sum is already zero. Right? In a certain sense, that's a gauge that you can propose. But this demanding this equation here, u transpose delta u equals zero, then this equation is automatically satisfied. Right? And if I assume this, then, and now we go toward the expressing this projection. If I assume this, then I can retrieve all of these factors from this equation. So I have equations for delta s, delta u and delta v from the delta y, which is then important for the projection. And let me explain maybe the first one, because it's simple. So if I take the delta y up here, I hit it from the left with your transpose, then this term cancels, your transpose delta u is zero. Here it cancels the u, it becomes an identity. And now I hit it from the right with v. Then what happens is that delta v transpose times v is equal to zero. So this one also cancels. So I'm just left with this term here. And the left one is the identity, the right one is the identity, so I have delta s. So there explains this formula down here. And the other ones can be explained in a similar way. Right? So I hit it from the right with v, one term cancels. And the other one, I can cancel with this one. So I can actually extract my variations in these three factors of the decomposition from the overall variation delta y. And the delta y is the thing I want to project, right? So in a certain sense here, I've expressed what my projection is onto the tangent space. And that's what I can do now. So coming back to the matrix differential equation, which looks like this, abstract matrix differential equation, I make this low-ranked ansatz, meaning I express y truncated at rank r as this singular value decomposition type factorization. And then I can insert this, I can project the right-hand side onto the tangent space. And what I end up with, very much in the spirit of what I've just shown, instead of variations delta s, delta u, delta v, I now get ordinary differential equations for s, u and v. So I get three differential equations for the three different factors. And this is why it's called dynamical low-rank. So I do not get the low-rank by doing singular value decompositions all the time. But the point is, and the nice thing about this, I get ordinary differential equations for all of the three factors. Especially what's interesting, I get ordinary differential equations for the bases. So the optimal bases in which I express the solution, they are collected in u and v, the orthogonal bases of the image space or whatever of my solution. They satisfy an evolution equation. So I get them in a certain sense automatically. The method selects basis functions by itself. All right. There are some ugly details to this, which I will just brush upon. But in principle, it's a very nice method. And now, coming a little bit to the motivation of this, if I look at now this system, if I now propagate these three factors of the decomposition instead of the full solution, what's my advantage? The advantage is that I only have a linear effort in the large discretization numbers. So think of my spatial discretization as n. Think of a 3D spatial discretization, a million degrees of freedom. Think of a velocity discretization in n, which might be in the hundreds or something like this. And typically, my whole quantity is of order m times n, right? So number of spatial cells times number of velocity. If I now look at these matrices, however, they are only low-ranked decomposition. So this is an n by r matrix. This is an n by r matrix. And this is an r square matrix. So the effort, if the rank is low enough, is only linear in each of the two. So linear in the number of spatial degrees of freedom plus linear in the number of velocity degrees. That's the big advantage. And even more important than computational effort is memory footprint these days. So the memory also reduces by a lot. It's also only linear in each of the unknowns. And by the way, you can do the same thing for space. If you're 2D in space, you can do another low-rank. And so you end up again with something that's linear, not something quadratic. That's the big advantage of this. All right. So now, ugly details. So it looks nice, but it turns out that this integrator is unstable. But there is actually a variant of this that is stable. So this has been investigated in the ODE community for a while. And I'm just going to write down this variant of the integrator for you so that you see how the algorithm looks like. So here's again our 1D plus 1D kinetic equation. Now coming back to that. And now I make just an abstract low-rank representation of f. And it's abstract on function spaces, mind you. So there is no discretization here. First, I argued with a discretization, but low-rank is actually, or dynamical low-rank is actually something that works on the continuous level in the first place. So I express my f with some basis functions in space. They are called x depend on t. So they change over time. My expansion coefficients, which I collect in this small matrix, are by r, s. And then I have basis functions for the velocity space. So an abstract low-rank decomposition of f on functional spaces. And then obviously I have orthogonal projections onto these spaces. So that's just a standard truncated Fourier series. And then, and again, here I brush a little bit over the details. I can again compute the projection onto the tangent space. And I can express it like this. It can be expressed in terms of these projections. And the method now, I made the font a little bit smaller so that it fits on one slide. The thing that makes the method stable, I'll try to explain, is that I do not propagate the factors u, s, and v separately, but I propagate u times s. Then I propagate s in a second step. And then I propagate s times v transpose. Don't ask me why there's a stable. This would take another talk, but that's a slight modification that ultimately makes this method stable. So in the first step, I have this decomposition. This k is actually u times s. And I can derive an ODE for this, which looks like this. Don't worry about the details. Then now I do need a qr decomposition or an svd decomposition because I need to retrieve the two factors again. So I take this k and I extract the x times s, which I can do using an svd. Then I propagate s and I get just an ordinary differential equation for s, which is r by r. I form what I call l, which is s times w times v in the previous notation. And I propagate the l, which is s times v. And then I do an svd again to extract the separate factors. So it's a slight modification of what I've just shown, but that's actually the algorithm that we put into practice here. Written down for a kinetic equation, using all of the notation angle integrals, space integrals are all in this discretization as well. All right. So that's how it looks like. Now what I do in order to put this into practice is I numerically represent these abstract basis functions. Now I can say I have my abstract dynamical low-rank basis functions, something in space, something in velocity, and now I discretize them somehow. For example, in space I take finite volumes, whatever you want to do. We know dg. And in velocity space, I also do something. I might use a Galerkin. I might use a collocation, something like this. And the thing that I need to do is I need to find efficient ways to compute these averages. So that's a tricky thing in order not to destroy the overall effort, but using quadrature rules, you can do this, actually. And then there's another point that I might explain maybe a year from now at a conference where we see each other, is that the hyperbolic discretization becomes a little bit tricky here, because I have this term in there. So this is a spatial integral of the x-basis functions and multiplied with the derivative of an x-basis function. So if you imagine that I express this in terms of finite volumes, I have something that is discontinuous, right? And I take the spatial derivative out of it. So I need to play with, you know, fluxes, and I need to think about stabilization here. We have something in the working, which is quite interesting, because you have to do something that you wouldn't think you would have to do, but I will be happy to report on this in a year from now. But you need to be careful about this. That's kind of the message. All right, so that works. And I'm going to show you a couple of numerical examples in a second. However, I also want to say that this type of method stays interesting, and I think for a lot of people in the community, because I hope I can convince you that it has potential, but it has a lot of drawbacks that need to be fixed. It has no positivity of the probability density function. I want that, of course, but I'm working with orthogonal bases, right? And if I have an orthogonal basis, it needs to be negative. Otherwise, it can't be orthogonal to the next basis function, right? So naturally, these basis functions are oscillatory. Think of them as, you know, sine functions or something like this in the Fourier series. So I have no guarantee that my kinetic solution will be positive, and I need to fix that. There's also conservation is not built into this scheme, right? Again, same reason. I'm using orthogonal basis here, and I have not built in conservation. I have not built in entropy in this thing, yeah? Again, so this is pure approximation theory based on low rank, so it's a total different philosophy of deriving these equations. So I don't have that, and I need to work in order to get this. And a final thing among many things that we like here in this community is asymptotic preserving, right? So it's Xi Jin's workshop here, so I need to say something about this. And that's also something that's not built in from the start, right? And we have many strategies to derive asymptotic preserving schemes, and I think it's interesting to combine them with them. And there are, I think, one or two works already out there of how to do that. So can I ask a question? Yeah. Well, I understand you have entropy because you don't have an entropy, but do you have like L2 stability, those kinds of things? Yes, yes. So there's a kind of L2 stability here, yes, because that's kind of, let's say, easy to achieve because it's an L2 framework, if you wish. So I'm working with orthogonality and scalar products. But actually, that I would say is also something new. And we use, so the thing that I just mentioned about the stabilization, it's actually an L2 stability analysis that gives us the stabilization that we need to do and CFL conditions and stuff like this. Thanks. All right. Well, so I've talked a lot about low rank. So I've not given you the motivation except for its low memory. So there is the question, why should we expect the solution to a kinetic equation to have low rank in the first place, right? And you can think about this for a second. There is a very simple and convincing answer to that. Before I show you this, and you're probably peeking at it already, that I think is something interesting to the analysts here in this room, because the rank of a solution is really hard to analyze. It's an object that you don't get from looking at a somal left space, right? The usual things we do to characterize solutions of partial differential equations, they are more related to sparsity, right? If the solution is in a high somal left space, it means my coefficients of a Fourier expansion go to zero fast, so I can set many of them to zero. But the notion of rank is something that's harder to kind of grasp with the analytical tools that are out there. So in that sense, it's also interesting. How can I see, on an analytical level, if a solution to an equation has low rank? But there is a striking argument in kinetic theory why we should expect the solution to have low rank, because we have asymptotic limits. So if you look at the equation that I've just written down, I put my usual epsilon small parameters in there, long time observation, small source, strong scattering. Then in the case of linear transport, as epsilon tends to zero, we know that our solution f tends towards something that is isotropic. It doesn't depend on velocity anymore. And that means if the solution loses the kinetic variable, it means it's rank one in that notion. So I can write it like this. I can write it with just a spatial basis function. I don't even need that coefficient anymore, and one in velocity space. So it's rank one. And the obvious hope is that if you're close to the asymptotic limit, your solution is of low rank. It might not be rank one anymore, but it's of low rank. But all of this is essentially an intuitive argument. The proofs would be interesting to see, actually, how can we actually characterize the rank. All right. So before I show you two numerical results, a little bit of a literature overview. This method actually is quite old. It was invented by Dirac for quantum dynamics in the 1930s. So 90 years ago. So I realized when I made these slides, I cannot just say 30 anymore, because we're close to 30 in the 2030. So it's almost 100 years old. It was basically discovered by mathematicians about 70 years later. So Lubisch, who was an ODE person, basically stripped it off the physics language and turned it into a numerical method. And that was done about 10, 15 years ago. And they have developed this method for ordinary differential equations. Then even more recently, like a couple of years ago, Lubisch and one of his postdocs realized that kinetic equations are a good application field for this. So they published something on Vlasov-Parsson in the first place. And then other people jumped in, including myself and others in the community. So starting with linear transport, there is an AP paper here, which I have in my list, which I should have listed here. And already people are starting to work on these issues that I've just mentioned, so conservation and all of these other things. All right, so just two numerical results to convince you that this method works. And it does actually also work far away from the asymptotic limit. So this is a test case that we like in linear transport. It's called checkerboard. And it mimics a nuclear reactor, if you wish. So what you have is a source of particles in the center. And then there are obstacles. And the obstacles look like on a checkerboard. So you see this here in the solution, you see where it's blue, there is actually obstacles. So it propagates through a field of obstacles, which look like a checkerboard. And it's a complicated test case. It usually, well, it's a torture test case, which shows the shortcomings of methods. That's why we like it. And at the bottom, you see kind of see a resolved solution, which however has artifacts, which are called re-effects. So you see these rays. Those are actually numerical artifacts, but you see kind of how the solution should look like. Time in the beginning. And that's your time evolution here. So it propagates out and it says all of these structures. And that was computed with a lot of degrees of freedom. So it's a 250 by 250 spatial grid. And it has 820 angular degrees of freedom. So quite a lot. You can do the math yourself, multiply these numbers. And then I show you here a rank 10 solution, and here a rank five solution. And you can also compute the reduction in memory and computational effort here. And you see that well, rank five doesn't seem to be enough, but rank 10 already captures the solution quite well. And it's a significant reduction, because I do have just 250 squared times 10 instead of 820. And the other evolution is just 10 times 820, not 250 squared times 820. So it's a significant reduction, and it looks quite good. Another torture test case that we like is the so-called line source test case. It's basically computing a Green's function, if you wish. So you have in the initial condition at time zero, you have a delta peak, which propagates outward into a medium. And it's again a torture test case. And you see something that's really highly resolved, which is now a spectral Galerkin method in Engel. They're called PN solutions. They usually show wave effects, wave patterns here. So it's like it looks more like dropping a stone into water than the actual analytical solution. But what we can do is now we can crank up the angular degrees of freedom. So this is P39, which is something like 1600 degrees of freedom in Engel. But at the same time, we can lower the rank. So this is a rank 210. But below it is a discretization that actually has 1600 degrees of freedom. And in a certain sense, dynamic or low rank picks the right, picks the right basis functions to actually have a solution, which is much closer to the original solution than this. And it has less degrees of freedom that it actually uses. All right. That brings me to the conclusion. I hope that I could convince you that that dynamic or low rank approximation is well suited for kinetic equations. I think there's a lot of work for everybody. And the thing that it seems to do, it even seems to work a little bit further away from the asymptotic limit. It gets more interesting if you're 3d plus 3d, because then you have hierarchical tensors and different hierarchical tensor formats. And there's a lot that can be done too. And the main selling point is it's a huge reduction in memory footprint and also computational cost. So that's the end of my talk. Thanks a lot for the attention. And I'm happy to take questions. So I have a question, if I may. So there is a class of method now, which people call tensor method. And it seems very connected to what you are doing. So are there subtleties that I am missing or is it in the class of a tensor method? Yeah. I don't know exactly what you mean by tensor method. So I would make a guess. I'm guessing you mean sparsity. Could that be the case? Like a sparse? Yes. But it's more than that. It's really to write, as you do, the solution as a sum of a simple product. So usually they put two. You put three, but usually they put two. And then they sum up. Yes. So if you do 3d plus 3d, that's what you would do. So then you would go to a tensor representation of your solution. Yes. But one thing I want to stress again. So there is a philosophical difference between dynamical low rank and by the concept of sparsity. So sparsity means I do a numerical discretization in some way. A huge tensor. So if I do a 3d plus 3d kinetic density, this is a six indexes. It's a six tensor. And then in the six tensor, I can do a hierarchical tensor decomposition. And there are many, many versions of this. And then I can demand sparsity. So then I can say I want most of my expansion coefficients to be zero. That's the notion of sparsity. And there are many methods to achieve this. There are versions of singular value decompositions and so on to do that. Dynamical low rank is a little bit different because it first of all works on the analytical level, right? I have a low rank decomposition of that. And I have actually evolution equations for my basis functions. So that's a slightly different thing. I hope that answers your question. Yes, absolutely. Thank you. Another question? Yeah, I have a question. Is there air estimating in terms of the rank? And how do you know what's the rank to stop? You know, you try like 100, 200. So that's a it's like our looking at shankation. Usually there's some framework for error analysis. Is that possible for this approximation? So, okay, so I'm just going to answer that question. So Lubisch, right? He's a really very good ODE guy, right? He wrote these books together with Haire and Wanner, you know, the numerics for ODE's bibles kind of. And there exists for ordinary differential equations, there exists one paper on an error analysis. And he told me it's the hardest proof he has ever done for ODE's. So I think for for partial differential equations, it's out of reach for now. But I think it's a very, very interesting field of research. That's my answer. And I think I'm too, I'm too stupid to do that error analysis. If he tells me that it's hard. But it's interesting. Yeah, yeah, I can send you that, I can send you that paper if you wish. I have a right, some of his papers in this direction. May I ask a question, Martin? Sure. Okay, what do you, what do your low rank extension become when you have boundary condition, for example, when you have boundary layer, for example? Oh yeah, thanks for that question. I should have put this onto my slides. That's difficult, obviously. Because what's really hard is to come up with this projection in together with boundary conditions. So we have actually no idea how to do that. And because you have to, well, we have an idea to do this. You can do this basically formally, right? But the question is, you would need to, for example, in space, I need to, let me go to the equations, right? So this is an equation that has a spatial discretization. And I have to come up with boundary conditions now for the, for the, for my spatial basis functions, right? Because, and that's a hard thing, actually. Yeah, so we can, we can basically take boundary conditions and project them. And that's, that's not, it's not a good idea. Yeah, so, so that's something that's hard. I mean, it works for, you know, vacuum boundary conditions, it works obviously for, for periodic boundary conditions. That's not a problem, but layers goes, goes on the list of things that need to be done for this. Okay.
 So thanks a lot for the invitation. I would have very much liked to be at CIRM, you know,  like everybody else here, but we'll have to do and there is hope that we will be able  to meet in person again in the future. And I guess all of us have realized what we're  missing. So let me talk about a new thing that I've been working on or my group has  been working on for a year with different collaborators who are listed here. And it's  I think potentially also bigger idea for numerical methods for kinetic equations especially.  So it's called dynamical low rank approximation. I will explain it using radiation transport  because that's what I've been doing, but it's a general method that I believe is of interest.  And I should say I'm not the only one or we're not the only ones who have realized that. So  let me explain what it is because it's nice mathematics. So for the moment, imagine a very  simple kinetic equation. And for the sake of simplicity, I'm going to assume that it's  just 1d plus 1d. So it's 1d in space. And it's 1d in velocity or in radiation transport.  It's an angle variable. So you have a probability density depending on time, one spatial variable  and one angle variable. X might be an R and the angle in radiation transport is usually  in the interval minus one one, but the details don't matter. And here I've written down  just just just the just the 1d kinetic equation. It's 1d linear transport. So you have advection  here this term. You have a scattering term here. And the fly is just the integral over  velocity space averaged. There is a total cross section here loss of particles absorption  scattering and there is a source here. But again, the form of the equation doesn't matter  for now. I've put everything on the right hand side because I'm going to consider this  to be an ordinary differential equation in a second in some functional space. Alright,  now assume that that you discretize this equation somehow or semi discretize it somehow  in in space and in velocity in these two one dimensional variables. And what I end up with  is is a quantity that depends on two indexes, right? The index I is the spatial discretization,  the index J is the angle of velocity discretization. Alright, so I I end up with a matrix. And  if I don't discretize in time, then this is a matrix differential equation that I get.  And in general, it can be any right hand side here for the next couple of slides. And it  can be nonlinear, if you wish. Yeah, so why is your discretized kinetic density? And this  is the right hand side you get. And obviously, this is this is something quite natural that  you get you had these two indexes in kinetic theory. And you can also imagine that these  matrices become quite large, right? So if you do a 3d plus 3d kinetic equation, you have  more indexes than this, but you can put all of this into into a matrix, right? And your  one direction is space and the other one is velocity. Alright, so so it's called dynamical  low rank. So what are low rank matrices that's your linear algebra one or two lecture? Well,  rank is most easily defined by a singular value decomposition. So that's what I've written  down here. So I'm speaking of a rank r approximation. So this y gets a little index r here. And  this looks like a singular value decomposition. So there is a matrix u times a matrix s times  a matrix v. And an SVD usually looks like this, right? So you have your object, which  should have been called y here. You have you have a matrix which looks like a towel, right?  So this is m by m. This is an m by r matrix, this s is r by r. And this is r by n. If the  original matrix was m by n and r is the rank. And in a singular value decomposition, this would be  a diagonal matrix. But for now, we allow these matrices also to be full. So that's if I can do  this now, then that defines a rank r matrix or a rank r approximation. If I take the full  singular value decomposition and I just cancel out some of the singular values, then I have a  rank r approximation of that matrix. Now, something very well known from your standard  linear algebra course. I should mention a related concept here, which is sparsity. So rank and low  rank is related to sparsity, but sparsity is a little bit different. Because in sparsity,  you would ask for a matrix, for example, this s matrix here to have a lot of zeros, right,  which is a different concept than asking for a low rank. So I will maybe discuss this in the end  and maybe there are questions about this. Because ultimately, the goal is to save some computation,  to save some memory. All right, so that's the definition of a rank r matrix. And the idea  is that we want to approximate this matrix from our discretization just by a lower rank matrix,  right, by a rank r approximation. And so that's what we do, right? So I view this as an ansatz,  right? I basically say I assume my solution is of low rank, which means it has this SVD-like  decomposition, and I insert it into the equation, into my abstract ordinary differential equation  for the matrix that I have. Obviously, this truncated or approximated low rank approximation  will not satisfy the equation anymore, right? If I insert it, I will end up with something  that's residual, something nonzero. And now comes basically the trick of this dynamical low rank  method. I do something that I also do in a Galerkin method. And I would explain it in such a way  that this is now a nonlinear Galerkin projection. So in Galerkin, you would say the residual I  demand is orthogonal to the ansatz space. That's the standard Galerkin thing that you do. Here,  we do something different. What we do is we take the right hand side and we project it onto the  tangent space of rank r matrices. Why do we do that? Well, this matrix doesn't satisfy the equation.  We need to make sure that we get an equation that it satisfies, and we need to make sure that we do  not run out of this space of low rank matrices, which is now not a linear space anymore, but it's  a nonlinear manifold. And so in this way, I make sure that if I start with a low rank matrix, if  I project the right hand side onto the tangent space, then it basically says the time derivative  of this thing is in the tangent space, then I remain for all times in the space of low rank  matrices. So that's why I need to project the equation. And again, this is a nonlinear object.  If I'm at one point in this manifold of low rank matrices, which is y r, then at this point,  there is a tangent space. And that's what I do for now. Now, the question is, what is the  tangent space? Can I characterize it? And can I express this projection in such a way? Because  ultimately, I want to design a numerical scheme. And the answer is yes, I can do that. And I want  to give you a short argument of how the tangent space looks like. It's not a proof, but it's the  essence of the proof. So if I have an SVD-like decomposition of a matrix, y is USV transpose,  then you can do a simple product rule, if you wish. So if this is a variation, delta y,  think of it as a derivative, then I have a product rule. I have the derivative of the  first one or the variation of the first one, the second one, and the third one here. So this would  be an element in my tangent space here. And then I can say something about these variations,  because I want u and v to be orthogonal. So that means u transpose u is the identity. If I differentiate  this, I get this equation here, namely that essentially my variation delta u times the u  matrix, the original u matrix has to be zero in a symmetric way here. So this equation.  And the same holds true for v, because it's also an orthogonal matrix. And for s, I just need that  the delta s is also an R by R matrix. All right. So then what I do, I make this a little bit simpler.  I satisfy these equations by demanding that each of these two terms in the sum is already zero.  Right? In a certain sense, that's a gauge that you can propose. But this demanding this equation  here, u transpose delta u equals zero, then this equation is automatically satisfied. Right?  And if I assume this, then, and now we go toward the expressing this projection. If I assume this,  then I can retrieve all of these factors from this equation. So I have equations for delta s,  delta u and delta v from the delta y, which is then important for the projection. And let me explain  maybe the first one, because it's simple. So if I take the delta y up here, I hit it from the left  with your transpose, then this term cancels, your transpose delta u is zero. Here it cancels the u,  it becomes an identity. And now I hit it from the right with v. Then what happens is that  delta v transpose times v is equal to zero. So this one also cancels. So I'm just left with this term  here. And the left one is the identity, the right one is the identity, so I have delta s. So there  explains this formula down here. And the other ones can be explained in a similar way. Right? So I hit  it from the right with v, one term cancels. And the other one, I can cancel with this one. So I can  actually extract my variations in these three factors of the decomposition  from the overall variation delta y. And the delta y is the thing I want to project, right? So in a  certain sense here, I've expressed what my projection is onto the tangent space. And that's what I can  do now. So coming back to the matrix differential equation, which looks like this, abstract matrix  differential equation, I make this low-ranked ansatz, meaning I express y truncated at rank r  as this singular value decomposition type factorization. And then I can insert this,  I can project the right-hand side onto the tangent space. And what I end up with, very much in the  spirit of what I've just shown, instead of variations delta s, delta u, delta v, I now get  ordinary differential equations for s, u and v. So I get three differential equations for the three  different factors. And this is why it's called dynamical low-rank. So I do not get the low-rank  by doing singular value decompositions all the time. But the point is, and the nice thing about  this, I get ordinary differential equations for all of the three factors. Especially what's  interesting, I get ordinary differential equations for the bases. So the optimal bases in which I  express the solution, they are collected in u and v, the orthogonal bases of the image space or  whatever of my solution. They satisfy an evolution equation. So I get them in a certain sense  automatically. The method selects basis functions by itself. All right. There are some ugly details  to this, which I will just brush upon. But in principle, it's a very nice method. And now,  coming a little bit to the motivation of this, if I look at now this system, if I now propagate  these three factors of the decomposition instead of the full solution, what's my advantage? The  advantage is that I only have a linear effort in the large discretization numbers. So think of my  spatial discretization as n. Think of a 3D spatial discretization, a million degrees of freedom.  Think of a velocity discretization in n, which might be in the hundreds or something like this.  And typically, my whole quantity is of order m times n, right? So number of spatial cells times  number of velocity. If I now look at these matrices, however, they are only low-ranked  decomposition. So this is an n by r matrix. This is an n by r matrix. And this is an r square matrix.  So the effort, if the rank is low enough, is only linear in each of the two. So linear in the number  of spatial degrees of freedom plus linear in the number of velocity degrees. That's the big advantage.  And even more important than computational effort is memory footprint these days. So the memory also  reduces by a lot. It's also only linear in each of the unknowns. And by the way, you can do the same  thing for space. If you're 2D in space, you can do another low-rank. And so you end up again with  something that's linear, not something quadratic. That's the big advantage of this. All right.  So now, ugly details. So it looks nice, but it turns out that this integrator is  unstable. But there is actually a variant of this that is stable. So this has been  investigated in the ODE community for a while. And I'm just going to write down this  variant of the integrator for you so that you see how the algorithm looks like.  So here's again our 1D plus 1D kinetic equation. Now coming back to that. And now I make just an  abstract low-rank representation of f. And it's abstract on function spaces, mind you. So there  is no discretization here. First, I argued with a discretization, but low-rank is actually, or  dynamical low-rank is actually something that works on the continuous level in the first place.  So I express my f with some basis functions in space. They are called x depend on t. So they  change over time. My expansion coefficients, which I collect in this small matrix, are by r, s.  And then I have basis functions for the velocity space. So an abstract low-rank decomposition  of f on functional spaces. And then obviously I have orthogonal projections onto these spaces.  So that's just a standard truncated Fourier series. And then, and again, here I brush a little bit  over the details. I can again compute the projection onto the tangent space. And I can express it  like this. It can be expressed in terms of these projections. And the method now, I made the font  a little bit smaller so that it fits on one slide. The thing that makes the method stable, I'll try  to explain, is that I do not propagate the factors u, s, and v separately, but I propagate u times s.  Then I propagate s in a second step. And then I propagate s times v transpose. Don't ask me why  there's a stable. This would take another talk, but that's a slight modification that ultimately  makes this method stable. So in the first step, I have this decomposition. This k is actually u times  s. And I can derive an ODE for this, which looks like this. Don't worry about the details. Then  now I do need a qr decomposition or an svd decomposition because I need to retrieve the  two factors again. So I take this k and I extract the x times s, which I can do using an svd.  Then I propagate s and I get just an ordinary differential equation for s, which is r by r.  I form what I call l, which is s times w times v in the previous notation. And I propagate  the l, which is s times v. And then I do an svd again to extract the separate factors. So it's  a slight modification of what I've just shown, but that's actually the algorithm that we put  into practice here. Written down for a kinetic equation, using all of the notation angle  integrals, space integrals are all in this discretization as well. All right. So that's  how it looks like. Now what I do in order to put this into practice is I numerically represent  these abstract basis functions. Now I can say I have my abstract dynamical low-rank basis functions,  something in space, something in velocity, and now I discretize them somehow. For example,  in space I take finite volumes, whatever you want to do. We know dg. And in velocity space,  I also do something. I might use a Galerkin. I might use a collocation, something like this.  And the thing that I need to do is I need to find efficient ways to compute these averages.  So that's a tricky thing in order not to destroy the overall effort, but using quadrature rules,  you can do this, actually. And then there's another point that I might explain maybe a year  from now at a conference where we see each other, is that the hyperbolic discretization becomes a  little bit tricky here, because I have this term in there. So this is a spatial integral of the  x-basis functions and multiplied with the derivative of an x-basis function. So if you  imagine that I express this in terms of finite volumes, I have something that is discontinuous,  right? And I take the spatial derivative out of it. So I need to play with, you know, fluxes,  and I need to think about stabilization here. We have something in the working, which is quite  interesting, because you have to do something that you wouldn't think you would have to do,  but I will be happy to report on this in a year from now. But you need to be careful about this.  That's kind of the message. All right, so that works. And I'm going to show you a couple of  numerical examples in a second. However, I also want to say that this type of method stays interesting,  and I think for a lot of people in the community, because I hope I can convince you that it has  potential, but it has a lot of drawbacks that need to be fixed. It has no positivity of the  probability density function. I want that, of course, but I'm working with orthogonal bases,  right? And if I have an orthogonal basis, it needs to be negative. Otherwise, it can't be  orthogonal to the next basis function, right? So naturally, these basis functions are oscillatory.  Think of them as, you know, sine functions or something like this in the Fourier series. So I  have no guarantee that my kinetic solution will be positive, and I need to fix that.  There's also conservation is not built into this scheme, right? Again, same reason. I'm using  orthogonal basis here, and I have not built in conservation. I have not built in entropy  in this thing, yeah? Again, so this is pure approximation theory based on low rank, so it's  a total different philosophy of deriving these equations. So I don't have that, and I need to  work in order to get this. And a final thing among many things that we like here in this community  is asymptotic preserving, right? So it's Xi Jin's workshop here, so I need to say something about  this. And that's also something that's not built in from the start, right? And we have many strategies  to derive asymptotic preserving schemes, and I think it's interesting to combine them with them.  And there are, I think, one or two works already out there of how to do that. So can I ask a question?  Yeah. Well, I understand you have entropy because you don't have an  entropy, but do you have like L2 stability, those kinds of things? Yes, yes. So there's a kind of  L2 stability here, yes, because that's kind of, let's say, easy to achieve because it's an L2  framework, if you wish. So I'm working with orthogonality and scalar products. But actually,  that I would say is also something new. And we use, so the thing that I just mentioned about  the stabilization, it's actually an L2 stability analysis that gives us the stabilization that  we need to do and CFL conditions and stuff like this. Thanks. All right. Well, so I've talked a  lot about low rank. So I've not given you the motivation except for its low memory. So there is  the question, why should we expect the solution to a kinetic equation to have low rank in the  first place, right? And you can think about this for a second. There is a very simple and convincing  answer to that. Before I show you this, and you're probably peeking at it already,  that I think is something interesting to the analysts here in this room, because the rank of  a solution is really hard to analyze. It's an object that you don't get from looking  at a somal left space, right? The usual things we do to characterize solutions of partial  differential equations, they are more related to sparsity, right? If the solution is in a high  somal left space, it means my coefficients of a Fourier expansion go to zero fast, so I can set  many of them to zero. But the notion of rank is something that's harder to kind of grasp with  the analytical tools that are out there. So in that sense, it's also interesting. How can I see, on  an analytical level, if a solution to an equation has low rank? But there is a striking argument  in kinetic theory why we should expect the solution to have low rank, because we have  asymptotic limits. So if you look at the equation that I've just written down, I put my usual epsilon  small parameters in there, long time observation, small source, strong scattering. Then in the case  of linear transport, as epsilon tends to zero, we know that our solution f tends towards something  that is isotropic. It doesn't depend on velocity anymore. And that means if the solution loses  the kinetic variable, it means it's rank one in that notion. So I can write it like this.  I can write it with just a spatial basis function. I don't even need that coefficient anymore, and  one in velocity space. So it's rank one. And the obvious hope is that if you're close to the  asymptotic limit, your solution is of low rank. It might not be rank one anymore, but it's of low rank.  But all of this is essentially an intuitive argument. The proofs would be interesting to see,  actually, how can we actually characterize the rank. All right. So before I show you two numerical  results, a little bit of a literature overview. This method actually is quite old. It was invented  by Dirac for quantum dynamics in the 1930s. So 90 years ago. So I realized when I made these  slides, I cannot just say 30 anymore, because we're close to 30 in the 2030. So it's almost 100 years  old. It was basically discovered by mathematicians about 70 years later. So Lubisch, who was an ODE  person, basically stripped it off the physics language and turned it into a numerical method.  And that was done about 10, 15 years ago. And they have developed this method  for ordinary differential equations. Then even more recently, like a couple of years ago,  Lubisch and one of his postdocs realized that kinetic equations are a good application field  for this. So they published something on Vlasov-Parsson in the first place. And then other people jumped  in, including myself and others in the community. So starting with linear transport, there is an  AP paper here, which I have in my list, which I should have listed here. And already people are  starting to work on these issues that I've just mentioned, so conservation and all of these other  things. All right, so just two numerical results to convince you that this method works. And it  does actually also work far away from the asymptotic limit. So this is a test case that we like in  linear transport. It's called checkerboard. And it mimics a nuclear reactor, if you wish. So what  you have is a source of particles in the center. And then there are obstacles. And the obstacles  look like on a checkerboard. So you see this here in the solution, you see where it's blue,  there is actually obstacles. So it propagates through a field of obstacles, which look like a  checkerboard. And it's a complicated test case. It usually, well, it's a torture test case,  which shows the shortcomings of methods. That's why we like it. And at the bottom, you see kind of  see a resolved solution, which however has artifacts, which are called re-effects. So you see these  rays. Those are actually numerical artifacts, but you see kind of how the solution should look like.  Time in the beginning. And that's your time evolution here. So it propagates out and it says  all of these structures. And that was computed with a lot of degrees of freedom. So it's a  250 by 250 spatial grid. And it has 820 angular degrees of freedom. So quite a lot. You can do  the math yourself, multiply these numbers. And then I show you here a rank 10 solution, and here  a rank five solution. And you can also compute the reduction in memory and computational effort here.  And you see that well, rank five doesn't seem to be enough, but rank 10 already captures the  solution quite well. And it's a significant reduction, because I do have just 250 squared  times 10 instead of 820. And the other evolution is just 10 times 820, not 250 squared times 820.  So it's a significant reduction, and it looks quite good. Another torture test case that we  like is the so-called line source test case. It's basically computing a Green's function, if you wish.  So you have in the initial condition at time zero, you have a delta peak, which propagates outward  into a medium. And it's again a torture test case. And you see something that's really highly resolved,  which is now a spectral Galerkin method in Engel. They're called PN solutions. They usually show  wave effects, wave patterns here. So it's like it looks more like dropping a stone into water  than the actual analytical solution. But what we can do is now we can crank up the angular degrees  of freedom. So this is P39, which is something like 1600 degrees of freedom in Engel. But at the  same time, we can lower the rank. So this is a rank 210. But below it is a discretization that  actually has 1600 degrees of freedom. And in a certain sense, dynamic or low rank picks the  right, picks the right basis functions to actually have a solution, which is much closer to the  original solution than this. And it has less degrees of freedom that it actually uses.  All right. That brings me to the conclusion. I hope that I could convince you that that  dynamic or low rank approximation is well suited for kinetic equations. I think there's a lot of  work for everybody. And the thing that it seems to do, it even seems to work a little bit further  away from the asymptotic limit. It gets more interesting if you're 3d plus 3d, because then  you have hierarchical tensors and different hierarchical tensor formats. And there's a lot  that can be done too. And the main selling point is it's a huge reduction in memory footprint  and also computational cost. So that's the end of my talk. Thanks a lot for the attention.  And I'm happy to take questions. So I have a question, if I may. So there is a class of method  now, which people call tensor method. And it seems very connected to what you are doing. So are  there subtleties that I am missing or is it in the class of a tensor method?  Yeah. I don't know exactly what you mean by tensor method. So I would make a guess. I'm guessing  you mean sparsity. Could that be the case? Like a sparse?  Yes. But it's more than that. It's really to write, as you do, the solution as a sum of  a simple product. So usually they put two. You put three, but usually they put two.  And then they sum up. Yes. So if you do 3d plus 3d, that's what you would do.  So then you would go to a tensor representation of your solution. Yes.  But one thing I want to stress again. So there is a philosophical difference between  dynamical low rank and by the concept of sparsity. So sparsity means I do a numerical  discretization in some way. A huge tensor. So if I do a 3d plus 3d kinetic density,  this is a six indexes. It's a six tensor. And then in the six tensor, I can do a hierarchical  tensor decomposition. And there are many, many versions of this. And then I can demand sparsity.  So then I can say I want most of my expansion coefficients to be zero. That's the notion  of sparsity. And there are many methods to achieve this. There are versions of singular value  decompositions and so on to do that. Dynamical low rank is a little bit different because it  first of all works on the analytical level, right? I have a low rank decomposition of that. And I have  actually evolution equations for my basis functions. So that's a slightly different  thing. I hope that answers your question. Yes, absolutely. Thank you.  Another question? Yeah, I have a question.  Is there air estimating in terms of the rank? And how do you know  what's the rank to stop? You know, you try like 100, 200. So that's a  it's like our looking at shankation. Usually there's some  framework for error analysis. Is that possible for this approximation?  So, okay, so I'm just going to answer that question. So Lubisch, right? He's a really  very good ODE guy, right? He wrote these books together with Haire and Wanner, you know, the  numerics for ODE's bibles kind of. And there exists for ordinary differential equations,  there exists one paper on an error analysis. And he told me it's the hardest proof he has ever done  for ODE's. So I think for for partial differential equations, it's out of reach for now. But I think  it's a very, very interesting field of research. That's my answer. And I think I'm too, I'm too  stupid to do that error analysis. If he tells me that it's hard. But it's interesting. Yeah,  yeah, I can send you that, I can send you that paper if you wish.  I have a right, some of his papers in this direction.  May I ask a question, Martin? Sure. Okay, what do you, what do your low rank extension become  when you have boundary condition, for example, when you have boundary layer, for example?  Oh yeah, thanks for that question. I should have put this onto my slides. That's difficult,  obviously. Because what's really hard is to come up with this projection in together with  boundary conditions. So we have actually no idea how to do that. And because you have to,  well, we have an idea to do this. You can do this basically formally, right? But the question is,  you would need to, for example, in space, I need to, let me go to the equations, right?  So this is an equation that has a spatial discretization. And I have to come up with  boundary conditions now for the, for the, for my spatial basis functions, right? Because,  and that's a hard thing, actually. Yeah, so we can, we can basically take boundary conditions  and project them. And that's, that's not, it's not a good idea. Yeah, so, so that's something  that's hard. I mean, it works for, you know, vacuum boundary conditions, it works obviously for,  for periodic boundary conditions. That's not a problem, but layers goes, goes on the list of  things that need to be done for this. Okay.
2021-03-25_Frank_0_47.png 2021-03-25_Frank_48_215.png 2021-03-25_Frank_216_1547.png 2021-03-25_Frank_1548_1643.png 2021-03-25_Frank_1644_5219.png 2021-03-25_Frank_5220_8663.png 2021-03-25_Frank_8664_12335.png 2021-03-25_Frank_12336_12383.png 2021-03-25_Frank_12384_12479.png 2021-03-25_Frank_12480_12659.png 2021-03-25_Frank_12660_17303.png 2021-03-25_Frank_17304_21431.png 2021-03-25_Frank_21432_22211.png 2021-03-25_Frank_22212_22391.png 2021-03-25_Frank_22392_25043.png 2021-03-25_Frank_25044_27743.png 2021-03-25_Frank_27744_34103.png 2021-03-25_Frank_34104_34343.png 2021-03-25_Frank_34344_34763.png 2021-03-25_Frank_34764_38075.png 2021-03-25_Frank_38076_40379.png 2021-03-25_Frank_40380_43355.png 2021-03-25_Frank_43356_45515.png 2021-03-25_Frank_45516_46631.png 2021-03-25_Frank_46632_46751.png 2021-03-25_Frank_46752_48059.png 2021-03-25_Frank_48060_48167.png 2021-03-25_Frank_48168_53039.png 2021-03-25_Frank_53040_53171.png 2021-03-25_Frank_53172_53231.png 2021-03-25_Frank_53232_53399.png 2021-03-25_Frank_53400_53447.png 2021-03-25_Frank_53448_53519.png 2021-03-25_Frank_53520_53615.png 2021-03-25_Frank_53616_53855.png 2021-03-25_Frank_53856_53915.png 2021-03-25_Frank_53916_54023.png 2021-03-25_Frank_54024_54071.png 2021-03-25_Frank_54072_54143.png 2021-03-25_Frank_54144_54227.png 2021-03-25_Frank_54228_54407.png 2021-03-25_Frank_54408_54527.png 2021-03-25_Frank_54528_54635.png 2021-03-25_Frank_54636_54695.png 2021-03-25_Frank_54696_54815.png 2021-03-25_Frank_54816_54899.png 2021-03-25_Frank_54900_54947.png 2021-03-25_Frank_54948_55091.png 2021-03-25_Frank_55092_55163.png 2021-03-25_Frank_55164_55211.png 2021-03-25_Frank_55212_55259.png 2021-03-25_Frank_55260_55326.png
So thanks a lot for the invitation. I would have very much liked to be at CIRM, you know, like everybody else here, but we'll have to do and there is hope that we will be able to meet in person again in the future. And I guess all of us have realized what we're missing. So let me talk about a new thing that I've been working on or my group has been working on for a year with different collaborators who are listed here. And it's I think potentially also bigger idea for numerical methods for kinetic equations especially. So it's called dynamical low rank approximation. I will explain it using radiation transport because that's what I've been doing, but it's a general method that I believe is of interest. And I should say I'm not the only one or we're not the only ones who have realized that. So let me explain what it is because it's nice mathematics. So for the moment, imagine a very simple kinetic equation. And for the sake of simplicity, I'm going to assume that it's just 1d plus 1d. So it's 1d in space. And it's 1d in velocity or in radiation transport. It's an angle variable. So you have a probability density depending on time, one spatial variable and one angle variable. X might be an R and the angle in radiation transport is usually in the interval minus one one, but the details don't matter. And here I've written down just just just the just the 1d kinetic equation. It's 1d linear transport. So you have advection here this term. You have a scattering term here. And the fly is just the integral over velocity space averaged. There is a total cross section here loss of particles absorption scattering and there is a source here. But again, the form of the equation doesn't matter for now. I've put everything on the right hand side because I'm going to consider this to be an ordinary differential equation in a second in some functional space. Alright, now assume that that you discretize this equation somehow or semi discretize it somehow in in space and in velocity in these two one dimensional variables. And what I end up with is is a quantity that depends on two indexes, right? The index I is the spatial discretization, the index J is the angle of velocity discretization. Alright, so I I end up with a matrix. And if I don't discretize in time, then this is a matrix differential equation that I get. And in general, it can be any right hand side here for the next couple of slides. And it can be nonlinear, if you wish. Yeah, so why is your discretized kinetic density? And this is the right hand side you get. And obviously, this is this is something quite natural that you get you had these two indexes in kinetic theory. And you can also imagine that these matrices become quite large, right? So if you do a 3d plus 3d kinetic equation, you have more indexes than this, but you can put all of this into into a matrix, right? And your one direction is space and the other one is velocity. Alright, so so it's called dynamical low rank. So what are low rank matrices that's your linear algebra one or two lecture? Well, rank is most easily defined by a singular value decomposition. So that's what I've written down here. So I'm speaking of a rank r approximation. So this y gets a little index r here. And this looks like a singular value decomposition. So there is a matrix u times a matrix s times a matrix v. And an SVD usually looks like this, right? So you have your object, which should have been called y here. You have you have a matrix which looks like a towel, right? So this is m by m. This is an m by r matrix, this s is r by r. And this is r by n. If the original matrix was m by n and r is the rank. And in a singular value decomposition, this would be a diagonal matrix. But for now, we allow these matrices also to be full. So that's if I can do this now, then that defines a rank r matrix or a rank r approximation. If I take the full singular value decomposition and I just cancel out some of the singular values, then I have a rank r approximation of that matrix. Now, something very well known from your standard linear algebra course. I should mention a related concept here, which is sparsity. So rank and low rank is related to sparsity, but sparsity is a little bit different. Because in sparsity, you would ask for a matrix, for example, this s matrix here to have a lot of zeros, right, which is a different concept than asking for a low rank. So I will maybe discuss this in the end and maybe there are questions about this. Because ultimately, the goal is to save some computation, to save some memory. All right, so that's the definition of a rank r matrix. And the idea is that we want to approximate this matrix from our discretization just by a lower rank matrix, right, by a rank r approximation. And so that's what we do, right? So I view this as an ansatz, right? I basically say I assume my solution is of low rank, which means it has this SVD-like decomposition, and I insert it into the equation, into my abstract ordinary differential equation for the matrix that I have. Obviously, this truncated or approximated low rank approximation will not satisfy the equation anymore, right? If I insert it, I will end up with something that's residual, something nonzero. And now comes basically the trick of this dynamical low rank method. I do something that I also do in a Galerkin method. And I would explain it in such a way that this is now a nonlinear Galerkin projection. So in Galerkin, you would say the residual I demand is orthogonal to the ansatz space. That's the standard Galerkin thing that you do. Here, we do something different. What we do is we take the right hand side and we project it onto the tangent space of rank r matrices. Why do we do that? Well, this matrix doesn't satisfy the equation. We need to make sure that we get an equation that it satisfies, and we need to make sure that we do not run out of this space of low rank matrices, which is now not a linear space anymore, but it's a nonlinear manifold. And so in this way, I make sure that if I start with a low rank matrix, if I project the right hand side onto the tangent space, then it basically says the time derivative of this thing is in the tangent space, then I remain for all times in the space of low rank matrices. So that's why I need to project the equation. And again, this is a nonlinear object. If I'm at one point in this manifold of low rank matrices, which is y r, then at this point, there is a tangent space. And that's what I do for now. Now, the question is, what is the tangent space? Can I characterize it? And can I express this projection in such a way? Because ultimately, I want to design a numerical scheme. And the answer is yes, I can do that. And I want to give you a short argument of how the tangent space looks like. It's not a proof, but it's the essence of the proof. So if I have an SVD-like decomposition of a matrix, y is USV transpose, then you can do a simple product rule, if you wish. So if this is a variation, delta y, think of it as a derivative, then I have a product rule. I have the derivative of the first one or the variation of the first one, the second one, and the third one here. So this would be an element in my tangent space here. And then I can say something about these variations, because I want u and v to be orthogonal. So that means u transpose u is the identity. If I differentiate this, I get this equation here, namely that essentially my variation delta u times the u matrix, the original u matrix has to be zero in a symmetric way here. So this equation. And the same holds true for v, because it's also an orthogonal matrix. And for s, I just need that the delta s is also an R by R matrix. All right. So then what I do, I make this a little bit simpler. I satisfy these equations by demanding that each of these two terms in the sum is already zero. Right? In a certain sense, that's a gauge that you can propose. But this demanding this equation here, u transpose delta u equals zero, then this equation is automatically satisfied. Right? And if I assume this, then, and now we go toward the expressing this projection. If I assume this, then I can retrieve all of these factors from this equation. So I have equations for delta s, delta u and delta v from the delta y, which is then important for the projection. And let me explain maybe the first one, because it's simple. So if I take the delta y up here, I hit it from the left with your transpose, then this term cancels, your transpose delta u is zero. Here it cancels the u, it becomes an identity. And now I hit it from the right with v. Then what happens is that delta v transpose times v is equal to zero. So this one also cancels. So I'm just left with this term here. And the left one is the identity, the right one is the identity, so I have delta s. So there explains this formula down here. And the other ones can be explained in a similar way. Right? So I hit it from the right with v, one term cancels. And the other one, I can cancel with this one. So I can actually extract my variations in these three factors of the decomposition from the overall variation delta y. And the delta y is the thing I want to project, right? So in a certain sense here, I've expressed what my projection is onto the tangent space. And that's what I can do now. So coming back to the matrix differential equation, which looks like this, abstract matrix differential equation, I make this low-ranked ansatz, meaning I express y truncated at rank r as this singular value decomposition type factorization. And then I can insert this, I can project the right-hand side onto the tangent space. And what I end up with, very much in the spirit of what I've just shown, instead of variations delta s, delta u, delta v, I now get ordinary differential equations for s, u and v. So I get three differential equations for the three different factors. And this is why it's called dynamical low-rank. So I do not get the low-rank by doing singular value decompositions all the time. But the point is, and the nice thing about this, I get ordinary differential equations for all of the three factors. Especially what's interesting, I get ordinary differential equations for the bases. So the optimal bases in which I express the solution, they are collected in u and v, the orthogonal bases of the image space or whatever of my solution. They satisfy an evolution equation. So I get them in a certain sense automatically. The method selects basis functions by itself. All right. There are some ugly details to this, which I will just brush upon. But in principle, it's a very nice method. And now, coming a little bit to the motivation of this, if I look at now this system, if I now propagate these three factors of the decomposition instead of the full solution, what's my advantage? The advantage is that I only have a linear effort in the large discretization numbers. So think of my spatial discretization as n. Think of a 3D spatial discretization, a million degrees of freedom. Think of a velocity discretization in n, which might be in the hundreds or something like this. And typically, my whole quantity is of order m times n, right? So number of spatial cells times number of velocity. If I now look at these matrices, however, they are only low-ranked decomposition. So this is an n by r matrix. This is an n by r matrix. And this is an r square matrix. So the effort, if the rank is low enough, is only linear in each of the two. So linear in the number of spatial degrees of freedom plus linear in the number of velocity degrees. That's the big advantage. And even more important than computational effort is memory footprint these days. So the memory also reduces by a lot. It's also only linear in each of the unknowns. And by the way, you can do the same thing for space. If you're 2D in space, you can do another low-rank. And so you end up again with something that's linear, not something quadratic. That's the big advantage of this. All right. So now, ugly details. So it looks nice, but it turns out that this integrator is unstable. But there is actually a variant of this that is stable. So this has been investigated in the ODE community for a while. And I'm just going to write down this variant of the integrator for you so that you see how the algorithm looks like. So here's again our 1D plus 1D kinetic equation. Now coming back to that. And now I make just an abstract low-rank representation of f. And it's abstract on function spaces, mind you. So there is no discretization here. First, I argued with a discretization, but low-rank is actually, or dynamical low-rank is actually something that works on the continuous level in the first place. So I express my f with some basis functions in space. They are called x depend on t. So they change over time. My expansion coefficients, which I collect in this small matrix, are by r, s. And then I have basis functions for the velocity space. So an abstract low-rank decomposition of f on functional spaces. And then obviously I have orthogonal projections onto these spaces. So that's just a standard truncated Fourier series. And then, and again, here I brush a little bit over the details. I can again compute the projection onto the tangent space. And I can express it like this. It can be expressed in terms of these projections. And the method now, I made the font a little bit smaller so that it fits on one slide. The thing that makes the method stable, I'll try to explain, is that I do not propagate the factors u, s, and v separately, but I propagate u times s. Then I propagate s in a second step. And then I propagate s times v transpose. Don't ask me why there's a stable. This would take another talk, but that's a slight modification that ultimately makes this method stable. So in the first step, I have this decomposition. This k is actually u times s. And I can derive an ODE for this, which looks like this. Don't worry about the details. Then now I do need a qr decomposition or an svd decomposition because I need to retrieve the two factors again. So I take this k and I extract the x times s, which I can do using an svd. Then I propagate s and I get just an ordinary differential equation for s, which is r by r. I form what I call l, which is s times w times v in the previous notation. And I propagate the l, which is s times v. And then I do an svd again to extract the separate factors. So it's a slight modification of what I've just shown, but that's actually the algorithm that we put into practice here. Written down for a kinetic equation, using all of the notation angle integrals, space integrals are all in this discretization as well. All right. So that's how it looks like. Now what I do in order to put this into practice is I numerically represent these abstract basis functions. Now I can say I have my abstract dynamical low-rank basis functions, something in space, something in velocity, and now I discretize them somehow. For example, in space I take finite volumes, whatever you want to do. We know dg. And in velocity space, I also do something. I might use a Galerkin. I might use a collocation, something like this. And the thing that I need to do is I need to find efficient ways to compute these averages. So that's a tricky thing in order not to destroy the overall effort, but using quadrature rules, you can do this, actually. And then there's another point that I might explain maybe a year from now at a conference where we see each other, is that the hyperbolic discretization becomes a little bit tricky here, because I have this term in there. So this is a spatial integral of the x-basis functions and multiplied with the derivative of an x-basis function. So if you imagine that I express this in terms of finite volumes, I have something that is discontinuous, right? And I take the spatial derivative out of it. So I need to play with, you know, fluxes, and I need to think about stabilization here. We have something in the working, which is quite interesting, because you have to do something that you wouldn't think you would have to do, but I will be happy to report on this in a year from now. But you need to be careful about this. That's kind of the message. All right, so that works. And I'm going to show you a couple of numerical examples in a second. However, I also want to say that this type of method stays interesting, and I think for a lot of people in the community, because I hope I can convince you that it has potential, but it has a lot of drawbacks that need to be fixed. It has no positivity of the probability density function. I want that, of course, but I'm working with orthogonal bases, right? And if I have an orthogonal basis, it needs to be negative. Otherwise, it can't be orthogonal to the next basis function, right? So naturally, these basis functions are oscillatory. Think of them as, you know, sine functions or something like this in the Fourier series. So I have no guarantee that my kinetic solution will be positive, and I need to fix that. There's also conservation is not built into this scheme, right? Again, same reason. I'm using orthogonal basis here, and I have not built in conservation. I have not built in entropy in this thing, yeah? Again, so this is pure approximation theory based on low rank, so it's a total different philosophy of deriving these equations. So I don't have that, and I need to work in order to get this. And a final thing among many things that we like here in this community is asymptotic preserving, right? So it's Xi Jin's workshop here, so I need to say something about this. And that's also something that's not built in from the start, right? And we have many strategies to derive asymptotic preserving schemes, and I think it's interesting to combine them with them. And there are, I think, one or two works already out there of how to do that. So can I ask a question? Yeah. Well, I understand you have entropy because you don't have an entropy, but do you have like L2 stability, those kinds of things? Yes, yes. So there's a kind of L2 stability here, yes, because that's kind of, let's say, easy to achieve because it's an L2 framework, if you wish. So I'm working with orthogonality and scalar products. But actually, that I would say is also something new. And we use, so the thing that I just mentioned about the stabilization, it's actually an L2 stability analysis that gives us the stabilization that we need to do and CFL conditions and stuff like this. Thanks. All right. Well, so I've talked a lot about low rank. So I've not given you the motivation except for its low memory. So there is the question, why should we expect the solution to a kinetic equation to have low rank in the first place, right? And you can think about this for a second. There is a very simple and convincing answer to that. Before I show you this, and you're probably peeking at it already, that I think is something interesting to the analysts here in this room, because the rank of a solution is really hard to analyze. It's an object that you don't get from looking at a somal left space, right? The usual things we do to characterize solutions of partial differential equations, they are more related to sparsity, right? If the solution is in a high somal left space, it means my coefficients of a Fourier expansion go to zero fast, so I can set many of them to zero. But the notion of rank is something that's harder to kind of grasp with the analytical tools that are out there. So in that sense, it's also interesting. How can I see, on an analytical level, if a solution to an equation has low rank? But there is a striking argument in kinetic theory why we should expect the solution to have low rank, because we have asymptotic limits. So if you look at the equation that I've just written down, I put my usual epsilon small parameters in there, long time observation, small source, strong scattering. Then in the case of linear transport, as epsilon tends to zero, we know that our solution f tends towards something that is isotropic. It doesn't depend on velocity anymore. And that means if the solution loses the kinetic variable, it means it's rank one in that notion. So I can write it like this. I can write it with just a spatial basis function. I don't even need that coefficient anymore, and one in velocity space. So it's rank one. And the obvious hope is that if you're close to the asymptotic limit, your solution is of low rank. It might not be rank one anymore, but it's of low rank. But all of this is essentially an intuitive argument. The proofs would be interesting to see, actually, how can we actually characterize the rank. All right. So before I show you two numerical results, a little bit of a literature overview. This method actually is quite old. It was invented by Dirac for quantum dynamics in the 1930s. So 90 years ago. So I realized when I made these slides, I cannot just say 30 anymore, because we're close to 30 in the 2030. So it's almost 100 years old. It was basically discovered by mathematicians about 70 years later. So Lubisch, who was an ODE person, basically stripped it off the physics language and turned it into a numerical method. And that was done about 10, 15 years ago. And they have developed this method for ordinary differential equations. Then even more recently, like a couple of years ago, Lubisch and one of his postdocs realized that kinetic equations are a good application field for this. So they published something on Vlasov-Parsson in the first place. And then other people jumped in, including myself and others in the community. So starting with linear transport, there is an AP paper here, which I have in my list, which I should have listed here. And already people are starting to work on these issues that I've just mentioned, so conservation and all of these other things. All right, so just two numerical results to convince you that this method works. And it does actually also work far away from the asymptotic limit. So this is a test case that we like in linear transport. It's called checkerboard. And it mimics a nuclear reactor, if you wish. So what you have is a source of particles in the center. And then there are obstacles. And the obstacles look like on a checkerboard. So you see this here in the solution, you see where it's blue, there is actually obstacles. So it propagates through a field of obstacles, which look like a checkerboard. And it's a complicated test case. It usually, well, it's a torture test case, which shows the shortcomings of methods. That's why we like it. And at the bottom, you see kind of see a resolved solution, which however has artifacts, which are called re-effects. So you see these rays. Those are actually numerical artifacts, but you see kind of how the solution should look like. Time in the beginning. And that's your time evolution here. So it propagates out and it says all of these structures. And that was computed with a lot of degrees of freedom. So it's a 250 by 250 spatial grid. And it has 820 angular degrees of freedom. So quite a lot. You can do the math yourself, multiply these numbers. And then I show you here a rank 10 solution, and here a rank five solution. And you can also compute the reduction in memory and computational effort here. And you see that well, rank five doesn't seem to be enough, but rank 10 already captures the solution quite well. And it's a significant reduction, because I do have just 250 squared times 10 instead of 820. And the other evolution is just 10 times 820, not 250 squared times 820. So it's a significant reduction, and it looks quite good. Another torture test case that we like is the so-called line source test case. It's basically computing a Green's function, if you wish. So you have in the initial condition at time zero, you have a delta peak, which propagates outward into a medium. And it's again a torture test case. And you see something that's really highly resolved, which is now a spectral Galerkin method in Engel. They're called PN solutions. They usually show wave effects, wave patterns here. So it's like it looks more like dropping a stone into water than the actual analytical solution. But what we can do is now we can crank up the angular degrees of freedom. So this is P39, which is something like 1600 degrees of freedom in Engel. But at the same time, we can lower the rank. So this is a rank 210. But below it is a discretization that actually has 1600 degrees of freedom. And in a certain sense, dynamic or low rank picks the right, picks the right basis functions to actually have a solution, which is much closer to the original solution than this. And it has less degrees of freedom that it actually uses. All right. That brings me to the conclusion. I hope that I could convince you that that dynamic or low rank approximation is well suited for kinetic equations. I think there's a lot of work for everybody. And the thing that it seems to do, it even seems to work a little bit further away from the asymptotic limit. It gets more interesting if you're 3d plus 3d, because then you have hierarchical tensors and different hierarchical tensor formats. And there's a lot that can be done too. And the main selling point is it's a huge reduction in memory footprint and also computational cost. So that's the end of my talk. Thanks a lot for the attention. And I'm happy to take questions. So I have a question, if I may. So there is a class of method now, which people call tensor method. And it seems very connected to what you are doing. So are there subtleties that I am missing or is it in the class of a tensor method? Yeah. I don't know exactly what you mean by tensor method. So I would make a guess. I'm guessing you mean sparsity. Could that be the case? Like a sparse? Yes. But it's more than that. It's really to write, as you do, the solution as a sum of a simple product. So usually they put two. You put three, but usually they put two. And then they sum up. Yes. So if you do 3d plus 3d, that's what you would do. So then you would go to a tensor representation of your solution. Yes. But one thing I want to stress again. So there is a philosophical difference between dynamical low rank and by the concept of sparsity. So sparsity means I do a numerical discretization in some way. A huge tensor. So if I do a 3d plus 3d kinetic density, this is a six indexes. It's a six tensor. And then in the six tensor, I can do a hierarchical tensor decomposition. And there are many, many versions of this. And then I can demand sparsity. So then I can say I want most of my expansion coefficients to be zero. That's the notion of sparsity. And there are many methods to achieve this. There are versions of singular value decompositions and so on to do that. Dynamical low rank is a little bit different because it first of all works on the analytical level, right? I have a low rank decomposition of that. And I have actually evolution equations for my basis functions. So that's a slightly different thing. I hope that answers your question. Yes, absolutely. Thank you. Another question? Yeah, I have a question. Is there air estimating in terms of the rank? And how do you know what's the rank to stop? You know, you try like 100, 200. So that's a it's like our looking at shankation. Usually there's some framework for error analysis. Is that possible for this approximation? So, okay, so I'm just going to answer that question. So Lubisch, right? He's a really very good ODE guy, right? He wrote these books together with Haire and Wanner, you know, the numerics for ODE's bibles kind of. And there exists for ordinary differential equations, there exists one paper on an error analysis. And he told me it's the hardest proof he has ever done for ODE's. So I think for for partial differential equations, it's out of reach for now. But I think it's a very, very interesting field of research. That's my answer. And I think I'm too, I'm too stupid to do that error analysis. If he tells me that it's hard. But it's interesting. Yeah, yeah, I can send you that, I can send you that paper if you wish. I have a right, some of his papers in this direction. May I ask a question, Martin? Sure. Okay, what do you, what do your low rank extension become when you have boundary condition, for example, when you have boundary layer, for example? Oh yeah, thanks for that question. I should have put this onto my slides. That's difficult, obviously. Because what's really hard is to come up with this projection in together with boundary conditions. So we have actually no idea how to do that. And because you have to, well, we have an idea to do this. You can do this basically formally, right? But the question is, you would need to, for example, in space, I need to, let me go to the equations, right? So this is an equation that has a spatial discretization. And I have to come up with boundary conditions now for the, for the, for my spatial basis functions, right? Because, and that's a hard thing, actually. Yeah, so we can, we can basically take boundary conditions and project them. And that's, that's not, it's not a good idea. Yeah, so, so that's something that's hard. I mean, it works for, you know, vacuum boundary conditions, it works obviously for, for periodic boundary conditions. That's not a problem, but layers goes, goes on the list of things that need to be done for this. Okay.
numerical methods for kinetic equations velocity or in radiation transport. It's an angle equations for all of the three factors. dynamical low-rank is actually something that decomposition of f on functional spaces. And then all of the notation angle integrals, space integrals are all in this discretization as well. need to be fixed. It has no positivity of the has low rank? But there is a striking argument in this method for ordinary differential equations. resolved solution, which however has artifacts, approximation is well suited for kinetic huge reduction in memory footprint and also computational cost. So that's the end of my talk.