 Hey guys, welcome back schizone episode 18 topic today is going to be more on solving linear systems in particular using LU decomposition or factorization and Currently my cats are going nuts and the neighbors are boring to the center of the earth So a lot of distractions, but don't let that distract you because it's a very important video first question why are we bothering with with these systems of linear equations and the reason is Because almost everything like 90% of things in science and engineering are Linear or can be assumed to be linear and when you have multiple things Mutually interacting or somehow connected in the network You ideally want to solve all of these equations At once and that would be a system of linear equations and that's a huge part of structural analysis Physics simulations circuit analysis statistics and so many more things So it's extremely important and that's why we put a lot of effort into implementing this in a nice quick elegant way and This reminded me a lot of the tutorial island Interaction between the player and the chef tutor in old school run scape Didn't we do this before just now like lab one was about Kramer's rule. Didn't we do this? Didn't we solve linear systems? Yes, but that's not really a good way to do it There are much better ways Including the one in today's video Just to jog our memories. What is what does it mean to solve a linear system or system of linear equations? Well, here's a word problem from first grade Theater production charges $21 for adult tickets 15 for students They sell 102 tickets for almost two grand. How much do they make, you know, how many each type? What did they sell the idea is you make two equations equation one being that of the revenue So 21 a plus 15 s equals the proceeds and then equation two is that of the number of tickets a Plus s equals one or two So a is number of adult tickets sold as the number of student tickets and that is that these equations Somehow contain information that when you even if you were to plot them find the intersection You would see that They're both satisfied at that crossing point and you were supposed to solve for a and s So in this case the way I was taught the teeth to do this is I would take equation two multiply it by 15 and subtract from equation one That would remove the s term and Then I would solve for a 67 I plugged that back into equation two to solve for s 35 and then you can check your work and plug numbers back in and see That it it's the correct answer Alternately you can form this as a matrix problem like this on the right hand side here So this big matrix two by two is the co-appreciate matrix So you have 21 15 1 1 as elements you have your unknown vector and then you have the right hand side known from the equations Ultimately the same approach is valid for these equip for these Matrix problems as well. We just refer to these things a little bit differently We call this kind of operation here an elementary row operation and we refer to this you know plug-and-chug as forward and backward Substitution we'll talk about those more in a few minutes That's kind of where LU Decomposition comes in and what this is ultimately is just taking advantage of a particularly easy type of system to solve And that is a triangular system. So what does that mean? It means a system that has zeros above or below the main diagonal And why is that so nice because this basically means you can Just solve the full equation all the equations just with Substitution because you already have one of the unknowns Basically, and you can just take that value and plug it into other other equations and solve for the other unknowns So there's very little thinking required once you have a triangular system And so the idea behind LU is basically you factor a into two pieces You're not cutting in half. You're not dividing it in half You're just reconstructing two other matrices that have the form of lower and upper Triangularity that when you multiply them together you get a again Why do you do this? Well, it's because if you take LU and plug it into the six equation here a x equals b a Being the vector of coefficient or measure the coefficients x being your unknown vector and b being your known vector you Basically can plug in your factored components L and U and then you can solve this instead of a complicated dense system like a originally was you have two triangular matrices L and U and you can break, you know UX into intermediate unknown vector y and solve ly is b solve for y Again, this is easy because you have a triangular matrix here. We'll talk about what that means in a second And then once you have solved for y y being u times x you can solve for x so The overall process is this step one you factor that system that a into an L and a U upper and lower triangular matrix And we're gonna use a do a little Algorithm we're gonna have once along the diagonal of L For a couple reasons we'll talk about that later step two would be to use forward substitution to solve for the intermediate Vector y then you use backward Substitution to solve for the unknown vector x. That's your answer to the problem So step two and three are very easy step one is a bit more complex, but even still it's not all that challenging So what does it mean to solve these triangular systems? Why did I say forward and backward? Substitution is so nice Well, here's an example So so you have this matrix 1 0 0 3 2 0 1 4 3 lower triangular Times an unknown vector x y z equals some known right-hand side Now x is already solved for you know row times column x is 5 Right, that's just be obvious The idea is you know x is 5 take that value and Remove all contributions of x to any other Equations in this system. So what you do is I can subtract off three x's from Row two both in terms of the coefficients on the left and On the right so I'm going to subtract off three times x and what is x equal to five? So I'm taking three off the coefficient and three times x or three times five off the right-hand side Similarly for equation three here at the bottom row take one off the x Coefficient and one times x off of the right-hand side multiplication the product You do that you basically put zeros along the entire first column outside of the diagonal element Or I should say below diagonal element. So yeah, you've done something productive very quickly Now y is trivial to solve for if you think about it because if you take this row Times this column y is obviously equal to six You have like well six now you can subtract off contributions of y to the other lower Equations so in this case we're going to subtract off four times y from equation three Both in terms of the coefficient as well as the right-hand side y was six So 45 minus four times six you now have this equation At this point now z is trivial to solve for because you can take this row times this column obviously z is seven and Now you're done. So you can see very very simple math Simple algorithm can be implemented to solve a triangular system both lower triangular and as you could imagine Upper triangular as well. So in summary triangular matrices are very good Now that's great, but how do we get these fabled triangular matrices? How can we go from a fully dense matrix like this to a lower and upper? Factor well It's pretty easy There are there are algorithms you could use but ultimately it's pretty straightforward approach You could do it by hand because it's it's kind of like kind of like Sudoku in a way basically We're gonna assume first off the diagonal elements of this lower triangular matrix in orange are ones. I Think that makes our solution unique. I'm not sure about that, but I'm pretty sure And then the elements above the diagonal are going to be zeros on the contrary for the lower for the upper triangular matrix the elements are unknown and Elements below that diagonal are zeros So right off the bat you can immediately find the entire first row of the unknown matrix But let's start simple here So if you think about it the first row first column of this dense matrix has to be the product of the first row of L The first column of you That means that the unknown question mark entry here has to be one, right? If you think about row times column that has to be a one now you can solve for This question mark entry here because you know this one element is the product of this row times this column And so we know that that question mark has to also be a one Now what can you do? Well if you want to solve for This question mark entry right here You can do so by multiplying, you know this row times this column the zero kills that question mark And basically, you know that this question mark has to be a four You've got that now you can solve for This question mark right here If you think about it this row times this column has to multiply to negative one that question mark must be a negative one similarly the last question mark in that row Must also be a one for a similar, you know process this row times this column has to give you a one that entry must be a one Okay, great. How about the middle question mark here? Well, if you take the middle row times middle column That's supposed to equal to therefore that question mark must be a three Okay, great. Now you have a two at the third row second column of Is this matrix? How do you get that? Well, that means that this question mark needs to be a two And now this last column of the you know you matrix. How can we get Get this question mark here? Well take this row times this column the zero kills that question mark so you know that this entry for this question mark needs to be a two and lastly you can solve for last question mark picking the Last row times the last column has to be one Therefore your last unknown needs to be a negative seven and if you were to multiply this matrix times this matrix You would get this matrix So that's kind of the process. It's a very simple process You can do it by hand very easily to construct the L you factorization for any matrix You always have enough information. Well, you may not always but usually had enough information assuming there's not a bunch of zeros in this To to construct this factorization and in fact if we know that the diagonal elements of L are always going to be one We can just assume we know that and then that case you only have these elements and these elements To encode Because you know that these are all zeros these are all zeros and these are all ones in that case You can fit that entire construction into the original memory space So you can basically replace a with the factorization of a and you can basically embed those Coefficients into the original memory space. That's the idea However, if you're not organized about this you might delete or overwrite things that you need this approach that I did May not have been in the ideal order to preserve elements that we used And so there are algorithms you can use you can look them up. Here's the one that I implemented in the code today And this basically means that if you use this approach and you kind of go in this sequence Diagonally, you will never overwrite an element that you need For a future operation and you can still safely use the original memory space of the a matrix to encode your L you factorization and if you want to take a look that's in the linear algebra Directory under LU the composition in the soy hub suppository now last thing to talk about is what about this pivoting? Nonsense. Well, the question is if you look at my note here There's an operation in this algorithm that involves dividing by diagonal elements And so if that element is zero or if the element is very close to zero Your solution might start to break down And so ideally you want to figure out a way to avoid the diagonal elements from being zero or close to zero and the overall idea is that we can actually adjust the rows the equations in our Matrix left-hand side and right-hand side such that the elements on the diagonal are not zero in fact we can pick the Equations with the maximum absolute value in that Component of the of the matrix and so yeah dividing by diagonal elements You know can be bad if that's nearly zero and so we can swap rows around on both sides of the equation Such that that row has the highest absolute value in that column to be located in that diagonal and so the idea is you basically pre multiply your Two sides of your equation so if your equation was a x equals b We're going to deconstruct that a into l u but we're also going to pre multiply both sides by p p is a pivoting or permutation matrix and It kind of has the form of an identity matrix But a little bit adjusted sometimes and so yeah, you'd have p l u times x equals p times b and So basically if you didn't want to permute or pivot the rows in your matrix You would leave these To be on the diagonal, but if you if you wanted to you could basically pick Which rows you want in which order by pre multiplying this p matrix of the form like this So this would have swapped the first and second rows of your equations and If you want to check I'm not going to go into it too much in this video if you're curious you can check the PLU decomposition file in the code base With that out of the way we'll talk about the actual code And I won't go into too much detail with the code because it's just all a bunch of random assembly stuff But I'll talk about how we're going to use these functions to solve important problems okay, so if you go to the Suppository we'll have these four different examples. Well four different directories The first one is matrix types that will show us like well, I'll just let's go into it to see example 18a Let's just run it So what's this printing out? so it's got this pie matrix so everything in here is pie three point one three four by four matrix and We're basically checking is this matrix upper triangular lower triangular or diagonal And it's it's checking somehow it's checking and it's saying nope not upper nope not lower and no diagonal How about this one? This is pies? Above and on diagonal. So yes, it's upper triangular, but no, it's not lower triangular or diagonal How about this matrix? Well not upper triangular, but it is lower triangular. It's not diagonal. How about this matrix? It's all three this matrix is upper triangular or triangular and diagonal So how does this work? Check out the code. Let's see what it's doing. Oops. I shrink this a little bit What do we got here? We have some includes The I guess three that matter the most are These three that say is upper triangular is lower triangular and is diagonal their functions that I guess take inputs pointing to that matrix the size of the matrix and Some tolerance for your zero Maybe you would that you would assume that point zero zero zero zero one is also zero So yeah, this would be how you would check if your matrix was upper lower and diagonal and So pretty simple. It's just calling those functions Printing out whether or not the matrix is any of those three things Great, how about example B? This is the step-by-step process for L you Decomposition so let's let's look at the code first So what are we including here? We're including this L you the composition File and what does it do? It just takes an a matrix and what it's going to do is it's going to well? I can actually pull it up math linear algebra L you Decomposition This Uses a pivot list do the algorithm to deconstruct the square RSI by RSI matrix that address RDI Fill in precision numbers into lower and upper triangular matrices stored together in the original memory space so it's what I explained in the slides and So if we take a look at the code again We have that including also we're we have some functions to Implement forward and backwards substitution. We have some functions to just kind of copy things around and make An identity matrix just for the sake of a step-by-step So how does this how does this work? So? Let's run it What's going on here? So this is our original a matrix that we're trying to solve It's trying to solve a x equals b. Here's a and here's b So the first thing we do is we call that L you Decomposition matrix and this is the decomposed L you in place so in this case these elements these elements and these elements are all part of the you matrix and Then this element and these two are part of the L matrix Along with an implied set of three ones on the diagonal So here's how that actually looks here's you you can see you takes the Six terms from this deconstructed a matrix and then L takes those three terms Plus the implied ones on the diagonal Okay, so the new problem is a x equals l u x equals b and we can use force position to find what u x is I refer to that value as y in the slides, but you can solve using Forward substitution for that y vector and now your new equation is your new problem is u times x equals y And we can use backward substitution to find x in that case and here's what it looks like So this is the solution to our problem the problem a x equals b which is right here is Solved by that so I'm going to snag this and let's open up octave to check our work here Let's plug in a and b And if we take the solution of that system like this our answer should be negative point seven point six Negative one was that what it was let's run it again negative point seven point six and negative one so yes our Algorithm worked the l u Facturization worked as well as the backward and fourth institution also worked We effectively solved for the unknown vector x, so that's great. That's example b Now example c this was a more compact approach because you don't need to get all the intermediate steps all the time it's a lot of You know work, so in this case all we're doing is we're calling l u solve Let's look at the code this basically is a wrapper for everything else It will do all the work, but you have to pass in everything you have to pass in your destination vector for your unknowns you have to also pass in the A matrix that you want to factor and then use to solve the equation as well as your right-hand side vector Would be at the location in rdx as well as the size of your system square system of dimension rcx And so how does this work? Well, honestly, I think this what this is doing is let me run it for you So here is it's trying to solve ax equals b. Here is a large 10 by 10 System with a bunch of random numbers in it. You can see here. They're all between looks like zero and one And here is your right-hand side vector. And so we're trying to solve ax equals b. Here's a here's b and Here is the solved x. This is you know, we've called l u solve and it has spat out x Now we can actually check our work not only an octave, but even with our own program because we can we can see if Hold on let me clear everything If ax equals b and we have x well then I can subtract I can multiply our Computed value for x times our a matrix and Subtract off b and I should get I should get zero right so we can do that here I did that ax minus b is pretty much zero So let's see how this worked in code so the first thing that happens is a bunch of printing which is nonsense and then it calls l u solve with the proper input so it has a Memory address for the a matrix memory address for the b matrix memory address destination for the x matrix and the size of the matrix He's you call l u solve and it now has placed the answer in this nation at Rdi and so we have some place in memory here here dot x. That's a ten by one Array of doubles so eight bytes each To fit that unknown vector and so yeah pretty simple it it solves the problem for us in one step Very good. Now. What about example D? This one is a little bit special because it involves the pivoting And so it's not going to matter so much for this example because all the numbers were kind of all the same order and the matrix was fully populated but You know if you had numbers here that were very very close to zero or zero on the diagonal your pivoting would Be actually very useful. It's still going to have an effect on this But it's not going to be noticeable in terms of your answer because it it didn't really do anything valuable But anyway, it's the same entries for a and for b as before and in this case We've used instead of l u solve. We have another function called pl u solve Which not only spits out an x matrix Which is the answer, but it also spits out the permutation matrix So it enables you in this case you can see it wanted the instead of the zeroth row to be an entry zero It wanted the seventh row to be in the seventh to be in the first position Because that gave it the best The highest absolute value at the zero zero Diagonal it was in the seventh seventh row can we can we tell that let's see Looking at the first column of the a matrix is the seventh row aka the eighth row This one the highest absolute value the answer is yes And so you can begin to see how that permutation has an effect It's basically pulling the row that has the highest absolute value at the You know the ith column into the ith diagonal and so we've dragged the seventh row of our a matrix up here Similarly, we would then grab the fourth row the ninth row the sixth row You can see this p vector kind of tells you which rows we've grabbed in which order And remember you have to apply that operation not only to the a matrix but also to the b matrix Because it's not just one hand of the equation But both sides of the equation that you have to you know swap around to keep things Working, but again, even still you can see that our result is that a x minus b is zero meaning our solution for x is correct Lastly, let's check out the code here to see what it what it was doing. How is this different? Well, it looks like the only difference is we change the call of l u solve to pl you solve Which is all the same thing except I believe you pass in an additional Entry here. Yeah r8. I got to change that was there right now r8 needs to can use to contain the eight location that you want to return and Use to compute the permutation array and so let's put that in here right now PL you solve That's going to be a Double no, it's not it's an it's an int It's a u int At r8 So with that out of the way We've implemented not only l u but also pl u decomposition and we know that it works So we have the ability to solve these dense systems of any size So I showed you examples with 10 by 10 matrices, but you can imagine this extends to 100 by 100 matrices, but you can tell that it would obviously be much more expensive Imitationally with that I'll end the video. I want to thank you guys for watching if you want to hang out We have a discord server last link in the description. I'll see you there