Gaussian Elimination From Scratch in Python: A Step-by-Step Guide

Key Takeaways

what is going on guys welcome back in this video today we're going to implement gsan elimination from scratch and python which is an algorithm for solving systems of linear equations and implementing it from scratch can help you understand it better so let us get right into [Music] it all right so we're going to implement goshan elimination from scratch in Python today and the main purpose behind this is to understand the algorithm better it's not because we need it we can just use the linear algebra module from numpy if we want a solution to a system of linear equations the main purpose is to better understand goshan elimination by coding It Out by implementing the individual steps and uh we're going to do that in this video today now just as a quick recap I'm not going to teach you goshan elimination in general but for those of you who maybe learned about this some time ago just a very brief reminder for what uh goshan elimination is I'm going to open up my paint here briefly uh the basic idea is we have a system of equations something like 3x + 2 y - Z is equal to 5 and then we have something like 2x + y - uh 3 Z = 2 or something and we take such a system maybe we have third one here to make it complete uh 4x - 2 Z = 1 or something like this and we take this and turn it into a matrix so we have 3 2 -1 5 into an augmented Matrix here 2 y -3 or sorry not y uh 2 1 -3 2 4 0 -2 1 and the gsan elimination now basically is we can swap rows we can add uh equations and we can multiply with scalers and the goal is to get to the so-called uh reduced row Echelon form where we have once here and we have zeros here and then we can just uh substitute back the solutions and we get uh or we can do the substitution and then we can get uh the results for X Y and Z that's the very basic idea I'm going to maybe sketch it out again when we get to the uh individual sections but we're going to implement this from scratch in Python now so what we're going to use for this is we're going to use numpy but we're not going to use linear algebra from numpy so we're just going to use numpy to work with array not in order to solve the uh or to to do the algorithm for us so I'm going to open up a new file now let's go to my current directory let's remove everything that's in here and then let's create a main py file uh which we're going to just delete the swap file here so I'm going to zoom in and what I want to do first is I want to have a numpy based a linear algebra numpy based function that can give me the correct results so I can can double check if I'm doing everything right so I'm going to say import numpy SNP and then we're going to define the function solve system of equations and we're going to get a and b as input here a being the Matrix representing the coefficients and B being the target vector or results vector and we're going to just say now try to get a solution by saying numpy linear algebra solve a and b and then we can just return the solution and if we get the exception that is the numpy linear algebra Lin Al error so linear algebra error if we get that we return system has no unique solution so that is going to be our method that we use or our function that we use to double check if something is is uh done correctly or not we're going to use a simple example here first I'm going to Define an array and this array now is going to uh contain our individual row so we're going to say here 31 -1 that is equal to saying 3 * X - 1 * y - 1 * Z is equal to whatever is in B at that position uh then I'm going to say come on uh then I'm going to say one 1 0 I need to do something about this autoc completion and then we have two 0 -3 all right so this is our coefficient Matrix and then our results Vector is going to be NP array and it's going to just be uh 5 02 so these are our equations 3 or 3x - y - Z is equal = to 5 x + y is = 0 and uh x - 3 Z is = 2 and we're going to try to solve this system of equations now now I can do this right away with my function here so just a and b will give me the correct answer and this is what we're going to compare our solution with so 1.3 1.3 0.2 all right so let's solve this with Gan elimination implemented from scratch and for this I'm going to go ahead and say Gan elimination AB actually gsh and elimination is only the process of getting it into the reduced row Echelon form we can then also do the substitution which is I think not part of the gsh and elimination I might be wrong about this but we're going to do both we're going to basically solve the system um so what we're going to do first is we're going to create the augmented Matrix so I'm going to say n is equal first of all to the length of a to get the dimensions and then I'm going to say augmented is equal to NP and then we're going to going to use the column stack function um to just concatenate The Columns of A and B so we get this augmented Matrix with this uh line in between even though of course there is no line this is just some uh notation for humans uh and then we're going to turn this into a float Matrix so this is our augmented Matrix and now we're going to just uh perform the basic operation so what we want to do is we want to take the rows that have the highest pivot element and we want to swap the rows if necessary so that the row with the highest value for that pivot uh column that we're currently uh looking at is at the current position so maybe this was not the best way to phrase it let me show you visually what I mean by this um if we have some I'm not going to choose the one that we have now but if I have some Matrix that has maybe a value one two three in here and then also the other values here whatever they are uh and then I have maybe some three two two and then maybe I have uh 51 one then the goal is to get the highest pivot element the highest value here at this position so this is the pivot element and we want to have the highest value here so we would swap these two rows then for the second column uh for the second yeah for the second column we want to have the highest pivot element in this position so if this was for example uh three we would swap it so that the three is a disposition we want to have the highest value uh always at the pivot position so to say so we're going to do it the following way we're going to say for I in range n which is the dimensions of our MX I going to say that the pivot row by default is going to be just the first row and then we're going to say the maximum value at the position that we're interested in right now is just going to be that specific uh the value that we have so it's going to be the absolute value but by the way maximum means absolute maximum value so maximum absolute value we're not interested in positive values only want to have the highest absolute number there so maximum value is equal to Absolute value of augmented and then we're going to take the diagonal element in the first iteration it's going to be the first one then the second one and so on it's always the diagonal um and then we look for the highest value so we say four K in range I + 1 because I don't need to consider myself I don't need to consider the first row um up until n and we're going to say if the absolute value of that of the element in this row so augmented K still I because we're looking for the First Column then the second column then the third column if that is larger than the max value we're going to um say that the max value is equal to that augmented Matrix Ki value so that that would be if I use the example here again um it would be this value here so I would say this is augmented I I so 1 one or or actually 0 0 in the first iteration and then we would see okay this one is higher so now that would be my choice but then this one is even higher so that would be my choice uh so I basically keep track of the highest value which at the end would be five and I say that the pivot row is going to be K every time I find a new value so that is now my Pivot row and once I'm done with all the rows so I definitely found the best one what I'm going to do is I'm going to say if the pivot row did not remain the same so if it's not equal to I I'm just going to swap the row so I'm going to say augmented I um is going to be equal or augmented I augmented uh pivot row whatever it is it's just going to be swapped so we're going to say augmented pivot row and I'm going to use a copy here and augment mented I copy so we're just swapping The Columns so since this is the best uh since this is the highest value for the First Column this is my Pivot rle so I perform this swap so these two are now exchanged including of course these values so that is that and uh we're then going to check if this value it's the best one so we don't have to check for all the values but if this best value if this maximum value is close to zero uh I'm going to say that the system has no solution because that is um then basically uh problematic because if it's zero or very close to zero we cannot work with this so I'm going to say if my uh if the absolute value of augmented II after the swap so of my Pivot element which is the optimal one the largest one if that is below a certain threshold let's say I uh 1 E uh -10 which is very small I'm going to return system has no unique solution all right and otherwise I'm going to just continue I'm going to say that um that the pivot is equal to augmented I I so this is now in this case we would swap the rows and then this one would be here so I would have maybe actually let's do this to not keep this too abstract so I'm going to delete 511 and then here we had 133 I think so 51 1 1 3 3 and then we would have 15 and here we would delete the one so that would p a swap and then that would be my Pivot element so I would say augmented uh I I is the pivot and then I would say I want to divide this row by the pivot El so augmented I I shall become whatever it was before or actually the whole row sorry not just the pivot element uh the row should be whatever it was before divided by five so basically we're taking all these values here um and replacing them so this becomes one this becomes 1 over five this becomes 1 over five this becomes three and so on so we do this to get the reduced row Echelon form where we have the ones at the diagonals so we want to have one one one here so this why we do that uh we do this Division and then we do the following we eliminate the pivot element uh from from the rows below so remember row Echelon form reduced row Echelon form means 1 one one but it also means that uh these here have to be zero so we need to eliminate the three and the one so what we do is we say for J in range I + 1 to only consider the rows below what we do is we say give me the factor the factor is just the number that I have here so this is one so I only need to take this as many times as I have here so in order to eliminate uh this element here I just have to take this row minus that many times this row because I have a one here so in this case this row here minus three times this row here will result in this becoming a zero because 3 minus 3 is zero so factor is going to be the augmented ji the first element or the element in the same column but in all the other rows that's going to be the factor and then I just say augmented uh augmented J is going to be equal to augmented J minus the factor times augmented I which is of course valid we can add we can swap and we can multiply with scalers multiplication also means division adding also means subtracting so we can just do this and then we have the zeros below our pivot element and this of course is repeated now for all the different columns that thing that we just did here is repeated for is done once for this so after being done with that we're going to have a zero here and whatever we get uh when we take this minus three times that so we're going to have some other value here and here and here and uh then we're continuing with this so we're trying to find the pivot element uh or the row with the largest pivot element in this position it's going to be that so we swap these two rows and then we perform the same operations again and again until we have all the rols handled and that's actually it that is already now our Matrix in reduced row Echelon form so we already have the form done the only thing we need to do now to get the solution is we need to substitute back so we need to say um is this a proper indentation we need to say x is equal to np0 and size n uh and then we just say for I in range and think about it what do we need to do I'm going to use now a different example what do we need to do if I have a matrix that is already in reduced row Echelon form so if I have uh one one actually let's make this like this one one one then I have Zer here and then maybe I have something like two and three and two or let's maybe say one then I have some values here like 10 5 15 whatever it is this is not our example now but if I have this reduced row Echelon form what do I need to do well in order to get one of the variables I just need to get this one here and that's going to be equal to that so I know that Z is equal to 15 then I need to get this one here and subtract um subtract the value from uh from the other world so I need to take this but I also need to subtract one times this value so let me put this into code and then explain it again for I in range we're starting at uh n minus one and we're going until -1 which means we're going until zero but we need to say negative 1 because otherwise it's not included we go step size of negative 1 so we go from last row up to the previous rows and we say that XI is going to be equal to and then augmented and we want to have the last element of this row so we're taking uh whatever row we're at we're taking the laste element which is of course this always so that's the value but in order to get to that value in the last row we need to do nothing so by doing this for the last row we're already done uh but for the other rows we need to make sure that we actually subtract all the stuff that's still here so we need to get rid of all this so in order to say um what Y is we need to say okay Y is basically five but not exactly five because we need to subtract 1 * 15 from here to know what y actually is so we say for all the other rows below for J in range I + 1 until the end so for all the rows that come below the row I'm currently looking at uh we're going to say that we need to subtract from our value that we took um the value for that specific column times whatever that value was so so XJ in this case if we're already in this row we would say whatever augmented i j is so one times whatever XJ is and J is the one from before so J would be 15 so I would say it's 5 minus 15 uh and that's basically it once we're done with that we have solved the system of linear equations I hope I didn't make any mistakes but actually now if I run this I should get the same result I don't because I column stack uh I need to use a list right or Tuple okay we did some M we made some mistake why is that let me just double check what we or what I messed up here so that is correct that is correct that is correct pivot row pivot row I that should be fine that is fine uh oh that's a problem I hardcoded five why did I hardcode five of course it's the pivot whatever the pivot is there you go now it works so this is how you do Gan elimination now again the goal of this video is not to explain Gan elimination in general but to just briefly recap the code we augment we create the augmented matrix by uh combining the columns we go through all the rows we try to always pick the maximum value for the pivot element we basically the row so that the maximum element is always at uh the position of the pivot element and then we divide by that number to make it a one and then we subtract from all the rows below this row to make all the values below the pivot element zero and then what we do is we just um uh yeah basically that's it and then we we just uh substitute back in we take the last value which is Trivial because it's already written there and then we take the ones above but we subtract the values that are not zero in this row and then we get the solution that's how simple it is and to show you that this works also for larger matrices I have another example here that we can run I'm going to replace this one this is a larger example I can run this and you can see we get the same result so we implemented goshan elimination from scratch so that's it for today's video I hope you enjoyed it and hope you learned something if so let me know by hitting a like button and leaving a comment in the comment section down below and of course don't forget to subscribe to this Channel and hit hit the notification Bell to not miss a single future video for free other than that thank you much for watching see you in the next video and bye