Lecture 25: Saddle Points / Inf-sup condition | Video Lectures | Mathematical Methods for Engineers II | Mathematics

Download this transcript - PDF (English - US)

INTRODUCTION: The following content is provided by MIT OpenCourseWare under a Creative Commons license. Additional information about our license and MIT OpenCourseWare in general is available at OCW.MIT.edu.

PROFESSOR: OK. Should we go? Great. Well I hope you enjoyed Patriots Day, and now we'll go for the final weeks of the course. And I wanted start with just a thought about scientific computing, and the people who work in it, and the problems. So just sort of an orientation. And maybe my point is that in part of the subject, the situation started like just give me the matrix, tell me the equation. I don't want to know the physics behind it, just the form of the equation, and I'll discuss how it could be solved. And to contrast that with what you may often be meeting in your thesis even in project, the details, the physical consequences, the actual detail behavior of this solution is what you're looking for. So somehow there are those two parts of the world, and people who know understand the physical problem are on the second side. And people who maybe think just about matrices, like that matrix for example -- that saddle point matrix, are more on the first side. Like just give me the matrix and I'll figure out how to solve the linear system. If it's a large linear system, I'll invent preconditioners, study stability, whatever you have to do.

But today's lecture is kind of half and half, because that is the framework that we're in. But the problem comes from a real problem in fluids. And it's called the Stokes problem, not Navier-Stokes. So Stokes alone is dealing with this question when there is a diffusion term, and I've normalized the coefficient there, viscosity term, whatever we call it, to be 1. And there is the constraint of incompressability, constant pressure divergence p equals 0. I'm sorry not constant, divergence of p equals 0 is a constraint. First let's do this. So make some comment on the physical side of this equation. What we're missing that Navier-Stokes would put in is a convection term. A term of v dot grad v that takes account of the fact that there's a flow here, things are moving. We're not in a stationary situation. And that term has been dropped. So that's telling us then we are in the case of very, very slow flow, where that term could be neglected. And certainly neglecting that term makes our problem easier. It leaves us with this symmetric matrix. It doesn't make it totally easy as we'll see, but it's a lot easier than when the convection term -- which by the way is not symmetric and pushes outside this framework. If that's present then of course in computational fluid dynamics, that's the central problem. You have to deal with it.

So we're taking a model problem that doesn't have the convection, doesn't have that unsymmetric complexity. But actually it's become more and more important. It's more than a model problem now, because there are a lot of medical application. Say probably the flow of your blood pushed by the heart that's moving along at a decent rate. But flow into tiny capillaries is slow, and Stokes is a reasonable step. And then actually Stokes is also going to appear in Navier-Stokes calculations, because often when you have two difficult processes you do a kind of splitting method, and you drop one term for half a step. So you take for example, a Stokes step, and then you do what the fluid forces you to do. You move that process with the flow. OK. So where flowing slowly here. And the nice thing is that this fits our saddle point framework, that's why I'm calling this block matrix S.

And you see really it came from minimization. The minimize was the problem that leads to the Laplacion, that stands for the Laplace operator with a minus sign to make it positive definite. And then this constraint gets imposed by a Lagrange multiplier P. So P is the Lagrange multiplier here. Oh wrong. Divergence V is 0, incompressive. Now that's on the tape. It's now on the tape right. The vergence V is 0. OK. So you remember how our equations looked? I'm going to write the unknowns here, V and P, So that you see what everything is, and it equals f and 0. OK. So this is now S, so I chop off here. So there you see the equation in the saddle point form. So A is the gradient, and it breaks on P. And so P is a scalar. I better say if I put a little arrow here, V is a vector, it has to components say if we're in a SD problem; f has two components. P the Lagrange multiplier for this constraint, the constraint is just a scalar constraint the diverge is 0. You remember what divergence means; dv1 dx plus dv2 dy, so those are the two components of V, the vector. And its divergence is 0. And these are, well with the minus sign, are adjoints to each other. OK.

So what's new that makes this problem worth some study? I guess that up to now, we've taken problems like this and we've done elimination to get to a single equation. You remember that. So that single equation is A transpose CA. Let's see. We're multiplying that row by A transpose C and subtracting. Let me do the elimination under here. I multiply the first equation by A transpose C. Just give me the matrix mode for now. I multiply by A transpose C and I subtract. So the first equation is still there, C inverse v plus Ap is f. But the second equation now when I multiply by A transpose C, so I just put it up there, and subtract that produces a 0 there, it produced minus A transpose CAp equals minus A transpose Cf. It leads us to this crucial matrix. And in applications up to now that was the way to go. That was totally satisfactory. We eliminated v, we got an equation for p, we solved it, and then once we knew p we went back and found v. OK.

So why is this different? Why don't I want to get to A transpose CA, the stiffness matrix if we were in finite elements. Why do I want to stay with this two field system? Which in finite elements is going to be a mixed method of finite elements. I'm going to keep v and p in the problem for a simple reason, that C in this problem is awkward, C would be the inverse of the Laplacion. C is the inverse of the Laplacion d second by x square plus d second by the y square. That's what delta stands for. And that's an uncomfortable operator to work with, the inverse here. In matrix language when we discretize, it won't be sparse, right. The Laplacion is sparse, the Laplacion gives us the five point matrix, the matrix with five non-0's on every row - four on the diagonal and four minus 1's at matrix that we know, the K2D matrix. But this is its inverse, C is now its inverse. It's the inverse of that matrix K2D that we've study, and of course it's full. So taking this step produces a completely full matrix, and all of our efficient algorithms for large problems work for sparse matrices. So whereas this one is sparse. I mean this the Laplacion. So we have the five point matrix sitting here. Here we have a discreet gradient, and a discrete divergence. They're all derivatives, and you think of the derivatives as being replaced by differences.

All these matrices are sparse, but if we take that fatal step of eliminating to get here, the matrix is full. Because of this fact that the inverse of a sparse matrix tends to be a full matrix. And you don't want to work with that. So we don't want, that's my message, don't want A transpose A when this is not sparse. OK. So that's the point about this example, which in the notes on the website is the section called the saddle point Stokes problem. OK. All right. So now we have both a specific example of this Stokes equations, three equations really because this is two equations for v1 and v2 and this is the third equation, the constraint. OK. So we have both a specific example and a general case.

So we're now going to work with that matrix as a whole, we're not going to reduce it to A transpose C, although in principle we could. And the question is one of stability. So like stability is a question that didn't come up for our minimum principles when there were no constraints. And maybe I'll just put a little mention why. So this is the first time we've had to face the question of stability. What do I mean that S inverse is under control? I don't mean that it's easy to compute, the inverse of that block matrix. But I have to know that the matrix is invertible. And stability is always like one slight step more than just being invertible. It means it's invertible and you have some control of the inverse. You know, it's not just on the verge of disaster. So we want to know about the inverse of that matrix. That's essentially it, we want to know about the inverse of that matrix. And in the process, we really want to know about the inverse of this matrix. So we may not want to compute it, that was my point, but stability is telling us we have here a matrix that's not positive definite. That the key issue here. Our saddle point matrix is not positive definite. It has positive eigenvalues and negative eigenvalues. Let me just take some examples. OK.

So here's an S. So I have a positive thing there. Non-0 stuff there. There is a simplest example I could think up. That matrix is not positive definite. It's determinate is minus 1 actually. But I mean they're safely away from 0. it's not singular or close to singular either. So that two by two example no problem. Well I could study the invertibility of this particular matrix, that would be solving differential equations. In other words I could study the case RP and V under control if I have a bound on f. Well they are.

But the issue comes when I go to finite elements, when I look at a sub space. so have to try to make that clear. What I plan to do to solve a continuous problem is to discretize it. This is my continuous matrix, block matrix, symbolic of the differential equations. And it's OK, analysis would show that. But the question is if I approximate V by some combination of finite elements, let me call it Vh. If I approximate the pressure by some pressure element Ph, then I'm sort of projecting the problem down into these finite dimensional, finite elements sub spaces. And of course f will get projected down to fh. And the matrices, these operators, will turn into finite as far as matrices. And the question is that projection safe? So that really it. It's the stability. S inverse is under control, this is OK. But is the projection of S inverse under control? Is the projection of S. So if I project it, is the inverse of that under control? That's the question. What comes up in applications, this is a projection of the continuous problem, the differential equation, onto finite element spaces, finite element trial function. OK.

So we've spoken a little about finite elements in this course, and more in 18085, and you're probably have seen them elsewhere. You know that we can approximate velocities. Here would be a really important example, because it's an especially simple finite elements. If I take my region, suppose it happened to be that rectangle, I divide that rectangle into small squares, there's a small square. I'm using squares rather than triangles. And in each square I make some low degree approximation, low degree trial function, for the velocities and the pressures. And the question is what functions can I use? For example, one I would love to be able to use would be, so in each square I would like to be able to use piecewise constant pressures. Piecewise constant pressure, so ph equal constant in each element. Element is the square. So a different constant in each element. So that would give me my vector of phs. It would be one component for every little square. And now what about Vh? What trial function would I like to use for the velocities? Constant velocities wouldn't be good, that's further than I'm willing to go. I'll go piecewise linear in each square. Now piecewise linear in each triangle meant that it had to form a plus bx plus cy. Strictly speaking that's what I should mean by piecewise linear with a different coefficient abc in each piece but because we're in squares now, we've got four nodal values, and we want a fourth coefficient. This was appropriate on triangles, on a square we really include the bilinear term dxy. OK.

So we have four coefficienct which determine and can be determined from the four values of the node. You see that I'm describing the whole space. Now I cut the region into more squares, I've described a away of discretizing the differential difference equation, the weak form of continuous problem. When I plug in these functions would give me a matrix problem, a discrete problem, that will have this form. And the question is do I have any control? So it will give me this projected matrix now, finite dimensional instead of a continuous problem. And the question is Have I still got any control? Well the answer for this particular choice of pressure space and velocity space is no. Unfortunately the matrix S that appears instead, the C inverse h shall I say it would be like our K2D inverse. That was pretty close of what it would be. And the Ah and the Ah transpose h and the 0 that matrix has a famous null vector for this specific example.

I'm taking big jumps here, because I just want to say the main point. The main point here is that a null vector appears in this problem. S fails to be invertible, and that null vector is like finite element people know its name. It has a couple of names actually, maybe checkerboard. The null vector is a vector when the pressure is 1 minus 1,1 minus 1 in a checkerboard pattern over the whole region, and a corresponding velocity. Well the trouble is Ah times this P is 0. That's what kills you. It turns out that A times that particular pressure is 0. And if you lose control of A that way, you've lost control A transpose CA. You've lost everything.

So I guess coming back to my opening comments. Here's a case in which we have a physical problem, but it falls into a pattern, where the mathematics actually decide are we able to use the trial spaces that are most convenient. I say convenient rather than most accurate of course. If I wanted to improve the accuracy, I would go up in the degree, I would take pressures to be piecewise linear and velocities to be piecewise quadratic with more points. And I could check the stability or not, is the matrix invertible or not, and the answer might be better there. In fact the notes and any book on finite elements has a list. And I don't plan on making that list here. But it has a list of press spaces, pressure element and velocity elements, on which the first guy that I've mentioned would be like Q of 0, degree 0 and velocities Q1 degree 1, and the letter Q being used for quadrilateral If I divided it into triangles, I would use the letter p for polynomial in that case. And this is a fail. This one fails.

And then another list would be how about Q1, Q1? I don't know. I wouldn't bet much on Q1, Q1, and of course I can look at the notes to see what the answer. Or Q1, Q2 that's probably got a better chance. You know people are inventing finite elements all the time, and they're having to check the stability or not, the uniform invertibility or not of S. OK. I just want to say that this one is so attractive, because it's so simple that people don't give up on it. They want to go with this, this is a fail, but so stabilized. That one is so attractive and, you know, specifically that it's this checkerboard mode that's in the null space. So you just think, how am I going to get that mode out of the null? And the stabilizing method is you stick a matrix here for purely numerical reasons, and it'll be a matrix that's exactly has this checkerboard pattern, so that the checkerboard is not any longer in the null space. Maybe I call it minus some small number times some E matrix, it's often written E. And alpha you keep as small as you dare. But the effect is that it takes this checkerboard pressure out of the null space. You see the problem was that without that the columns of A were not independent. A times P could be 0. And if those columns are not independent you're lost. Let me put the 0 back. Why was that a bad thing? I want to now sort of connect you to the matrix theory.

If I have this particular pressure and AhP is 0, so that tells me that those columns of Ah are dependent. Then the columns of S are dependent, right. If I have something in the null space of Ah, then I just put a pressure, this pressure, in the null space of Ah, and I just take the v to be there 0, then these guys don't do anything. And that A times that bad pressure gives 0, so I have vp, a velocity pressure, that's in the null space. OK. So the way it gets fixed is to stick in some small stabilizing factor there. OK. All right now that's the specifics of this particular example. Well I didn't give all the specific. If I were really on the ball I would show in detail why that checkerboard appears and kills this choice. I would check each of these other choices, and beside each one I would put fail or success. And when it fails, I would try to identify the failure mode, so that if you wanted you could stabilize and get rid of that failure mode. OK.

What I'll do instead is to come to the general question here. So the general question is to deal with this matrix S. And I write it again C inverse A, A transpose 0. I want to look if it's invertible? Can I get some bound on its inverse? What do I have to know? And I'll just emphasize here, what's the real problem? Because we saw it right there. The real problem is if the column of A are dependent. In other words i.e. in other words, AP equals 0 has a non-0 solution. Then you're wiped out. I'll just repeat. If the columns of A are dependent, and then the columns of S are dependent, the matrix is singular, and you're finished. So somehow we have to get some muscle behind A. But remember A is a rectangular matrix. This C is m by m, this A that comes up in the constraint is m by n. It has, we hope, more rows than columns. But still the columns could be dependent, and in this situation they are. OK. So the question is we have to get some muscle behind A. And this is the famous condition so called in soup condition that you would find in Professor [? Batiste ?] finite element book. But maybe these notes kind of focus on the heart of the idea.

And again in some way it's conditioned on A really, and on the pressure spaces, and the velocity spaces that you choose. OK. So these you could really say are Ah, Ph, and Vh, but I'll skip writing h's on everything now. So I'm down in the discrete case, but I have a sequence of discrete problem which are growing as h gets smaller, and the question is do I have a fix beta that control this? And it turns out that this the requirement this in soup condition says that for every pressure there is a velocity city, every discrete pressure now, there is a discrete velocity such that this holds. So it's sort of forcing A to stay away from 0. This beta is a fix number, and this sort of is the right normalization that goes with it. So you see why v transpose A and P? You know if I look at the quadratic form, v transpose P transpose, you know, the quadratic form with the matrix is always multiply by a vector here multiplied by its transpose here, where is the A term? A will multiply P, and it'll be multiplied by v transpose. In the quadratic part, in the energy, is this v transpose AP. And that's what we can't let go to 0. Because If that part fails us, if we lose these and we have a 0 there, this part on the sub matrix can't save us. Right.

So before I say anything more about that. let me just emphasize what stability was in the symmetric case. So suppose I have a positive definite symmetric matrix. it's got somebody eigenvalues, lambda 1 to lambda n. OK. What you see on this board is going to explain why stability is automatic if I'm dealing with positive definite things. So suppose I project? Let me just project in the direct way. I'll throw away the last row and column. Suppose I throw away the last row and column. That part is now gone, so these are the eigenvalues of K. And of course being positive definite, they're bigger than 0. OK. So that was the whole matrix, but now what I'm going to do is I'm going to project on the first n minus 1 rows and columns. So I'm going to look at the sub matrix, K sub n minus 1. And it has some eigenvalues. Are they the same as those? No way. And you see I'm not in this situation, this was not positive definite. So I'm looking here at the positive definite case. This has some eigenvalues mu1 up to mu n-1. These will be the eigenvalues of the sub matrix, Kn minus 1. And linear algebra is just like the right place to look for some information on the connection between the mus and lambdas. They're certainly not equal.

Our first question is if I take a positive definite matrix with positive eigenvalues, I throw away the last row and column. Do I still have a positive definite? Is that sub matrix positive definite? And the answer is yes. The eigenvalues are all positive. And more than that, how about that smallest eigenvalue? Because the eigenvalue is measuring like how close are we to disaster, to singular. Can I say anything about that smallest eigenvalue? And the answer is yes. The smallest eigenvalue of Kn minus 1 is never below the smallest eigenvalue of the full matrix. So if the full matrix is OK, all the projections, all the sub spaces, all this sub matricies are even more OK. And up at the other end, at the top, the top eigenvalue in sub matrix never passes the top eigenvalue in the big matrix. So you see at both ends I have what I want, I have control. If I have some range of eigenvalues for the big matrix, then all its projections are inside that range. And therefore the condition number, which would be the ratio of this to this, is smaller than the condition number of the whole matrix. So the condition number of Kn minus 1 is less or equal the condition number of Kn. OK. And that's exactly why working with positive definite matrices is so comfortable. The stability is just automatic.

But here we're not working with positive definite matrices. The stability is not automatic, and we need this in soup condition. And I don't know whether to try. Maybe I'm a little hesitant to go through the steps that lead from there to there, because this lecture is already asking you to put a lot of things together without details. And the notes do the best job I could do in making those steps simple and clear. And this is what you actually check. You check this condition for any proposed finite element, and any different problem. It wouldn't have to be the Stokes problem. Other problems would give us a constraint matrix A. A sort of stiff matrix C inverse. And we could check the condition again. OK. So maybe I'm going to leave that statement for the notes. OK. So that's this in soup condition.

And of course, you see right away that if A had dependent columns, which is the big risk, if A had dependent columns that would mean that Ap is 0 for some p. And then there would be no v. No v could possibly be found that would save it. If A had dependent columns, that's bogey man that we're avoiding above all. If A had dependent columns, then there would be a solution to Ap equals 0, this left side would be 0, the right side wouldn't, and the in soup condition would fail. But the in soup condition tells us well what about if there are pressures, that where Ap get smaller and smaller as the mesh size changes. And the question is can I find a v and so on. So that's the full condition but not to lose sight of the fact that the central question, the key question is the independence or dependence of the columns of A. OK. Good. So let's see, maybe that's as much as I want to say about the Stokes example, and what it leads to in a finite element world. I guess I'm going to stay for another lecture or two with these minimizations. So maybe I'll just anticipate slightly, the start of the next lecture. Because I'm quite pleased about getting this next topic into the series of lectures.

Let me say what it is. So it's going to be on minimization of an Au minus b square. OK. Just as this lecture is seen, the possibility comes up at the columns of A not independent, that this problem is in trouble. So I want to regularize this problem. So this is a major application in which A, maybe A transpose is the easiest thing to say, A transpose a is nearly singular or just plain singular. I'll say or singular. How do we rescue the problem? Because numerically and theoretically it's a disaster as it stands. OK. It's ill posed. So let me introduce that word ill posed problem. And let me mention, since I'm just using these last few minutes as a kind of getting the big picture, where do ill posed problems come from, where does these things come from? They come typically from what are called inverse problems.

So the application is to inverse problems. What's an inverse problem? An inverse problem for a differential equation is you know some inputs and outputs, you know some solutions, but you don't know the equation. You don't know the coefficient. So you're trying to determine what are the material properties by looking to see what the outputs are for experimental inputs. And very frequently you don't have enough inputs to really determine what these -- so it's inverse problem to find the coefficients of the equation given some information on its solutions. some measurements. And of course they're probably noising measurements. I mean this is an enormous source of applications growing an importance, and how do we deal with it. Well your regularize it. So this is typical of this type of problem. You add on some multiples of something that's safe, like for example u square. So the results is here that this A transpose a that comes from this problem there's an alpha times the identity coming from this problem. I mean it's an ordinary least square problem, and it leads us to the matrix A transpose a, comes from the first part. But I'm just talking about the standard normal equation. It will produce an A transpose a from here. But from here it will produce an alpha times the identity.

And of course just the way this minus E worked up there, here I'm doing it even more directly, I'm pushing the problem away from 0. This alpha, this parameter and the choice of that parameter, is the million dollar question in this whole subject. But it's purpose is to regularize the problem, to take his nearly singular matrix and push it up. See this matrix, A transpose a, never has negative eigenvalues, the problem is 0 or very small, small positive. So by adding on alpha I, I push those away from 0. Over there we had the problem of eigenvalues being positive or negative, and by pushing in one direction we might just bring a negative eigenvalue right on to 0. That's not our problem in this upcoming lecture. In this upcoming lecture we know that nothing is negative, but we need to make it properly positive. OK. So that's what's coming.

Can I just say the application that I want to describe is the first one I've given in life sciences, and I hope not the last. I don't know if any of you can suggest, and if you can an email will be very much appreciated. Applications of this whole framework in like sciences. So you'll see a particular application of this ill posed problem in trying to determine gene structure from expressions.

OK. So that's what's coming Friday. And I realize here I'm talking about material that probably is not going to go in your projects, because of the timing. I'm thinking that you're projects are more likely to build on what you had inn project one, and the basic material about solving large systems and basic minimizations. OK. Thank you. See you Friday.