 So what are the basic operations you can do on matrices? Well, many of them you are very, very familiar with. And rewriting them in C++ is quite trivial. So I'll just flash the code, and then we'll move on if it's clear enough. So if I want to initialize a square matrix to identity, what's the identity matrix? 1 on the diagonals, 0 away from the diagonals. So I say for int Rx equal to 0 through rows, int Cx equals 0 through columns. Demat Rx Cx is equal to, if Rx equal to Cx, then 1, otherwise 0. So this is quite trivial. Everyone happy with this example? Creating an identity matrix? How do I initialize an upper triangular matrix with 1s? So my matrix is a square matrix. Rows is equal to columns. Now again, I set up the outer row loop. I set up an inner column loop. And inside, I have to now check if Rx is less than equal to Cx. So remember the diagonal and things above the diagonal have to be 1s. Things strictly below the diagonal have to be 0s. This is Rx. That is Cx. In this part, Cx is larger than Rx. On the diagonal, Cx is equal to Rx. So the initializer will look like Rx less than equal to Cx, then 1, otherwise 0. So that's how you initialize an upper triangular matrix with 1s. How do you transpose a square matrix? So the important thing is not to double transpose. Transposition happens across the diagonal from the two sides of the diagonal. And I only need to consider each pair once. So I set up the row loop between 0 and rows. But now the column loop doesn't need to touch from 0 through rows. That is below the diagonal. It only has to start from Rx plus 1. That keeps it above the diagonal. And once I get a pair Rx, Cx that I want to flip, Rx, Cx with Cx, Rx, then the rest is very trivial. So I say float temp equal to this. F mat Rx, Cx is F mat Cx, Rx, and F mat Cx, Rx is temp. So here is a pictorial demo. So no need to transpose diagonal or below redundantly, just do it once. And swap Rx, Cx with Cx, Rx. So here is the picture. I start off with the first matrix. I don't look at the diagonal 0, 0, no need. I look at 0, 1, and 1, 0. I swap those. Then I look at 0, 2, and 2, 0. I swap those, 2, 7, 3 like this. Then I look at, I don't need to look at 5. I only need to look at this last number, 1, 2, and 2, 1. 6 and 8, and I transpose those. So this 2D matrix is the transpose of that 2D matrix, 4, 7, 8, 4, 7, 8, 2, 3, 6, 2, 3, 6. So that's why you only need to do Cx equal to Rx plus 1. You shouldn't go from Cx equal to 0. That's a bug. Very simple examples. Now perhaps the first somewhat interesting example is matrix vector multiplication. So I have this matrix called A, rows times columns, not square in general. I have an input vector x with calls columns. And I want to output y with rows columns. So the picture looks like this. So I can have a fact matrix A like that. I can have a tall x with calls columns. And the output is y, which is shorter. So those are the dimensions. And how do I fill it up? I run across the columns of A. I run down the rows of x, multiply, and add. So the important thing is to look at it in pictures. I want to fill up y0 and y1 starting with A and x, which will appear gradually. So first, A00 and x0 are multiplied, stored in y0. Then they disappear. A01, x1 appear. Those are multiplied and added on to go inside y0. Finally, A02 and x2 are added together to go into y0. Then the second row appears. I start scanning x0 again. A10 x0 plus A11 x1 plus A12 x2. So that's the value. I sweep across the column of A, row by row, and I fill up the matrix. So coding it is very easy. So for example, in the beginning, rx is equal to 0. I initialize y0 to be equal to 0. And then I sweep across the columns one by one. And I accumulate it into y. So yrx plus equal to arxcx plus xcx. So how many people are comfortable with matrix vector multiplication? Show off hands, please. In that corner? So very easy stuff. And what can you do with matrix vector multiplications? Well, we can do a lot. So many of you will be aware that, suppose I take a 2 by 2 matrix, like 2, 1, 1, 3. That is nothing but a transformation of an input 2D space to an output 2D space. For example, I can take the vector x, which I'll call the preimage. If I multiply x0, x1 by 2, 1, 1, 3, I will get an output vector y, which we call the image. So the preimage, when multiplied by A, results in the image, which is y. And in general, suppose the brown vector here is x, the blue vector will be y. Now let's do the following thought experiment. So suppose I have the 2D Cartesian space. I consider this unit circle on it. And I take, say, this vector. This is my x. Suppose I multiply it by A to find y. Maybe my y looks like this. Suppose I plot in this dotted line, which connects the tip of x, the preimage, to the tip of y, which is the image. And suppose I do this while rotating x all around. We'll see where the tip y goes. In particular, we will take this matrix 2, 1, 1, 3 and see where that takes us. So this I haven't written as C++ code because I don't have graphics wired up well. I've written it as some messy looking SyLab code, which you don't need to even see. But let me execute this code. So that line connects the tip of x to the tip of y. So the tip of x describes a circle. I've designed that. What is the tip of y describing? Looks like an ellipse. It is an ellipse. So the tip of x is inside. The tip of y is the outer ellipse. And the array involved is 2, 1, 1, 3. Now let me put back a few other things. For example, let me put back x and y itself. The picture will get a little busier. But now that I've started with the sparse picture, you'll understand what's going on. So in a different color, I'm going to show the other vectors. So now I'm actually connecting the red rays to the y's, to the edge of the ellipse. So now we see the tips connecting from the tip of x to the tip of y, and also origin to y. So in summary, you can call this the eye of the storm or something. That's the picture. So things are happening like this. This point was, say, going to that point. So this point was going to that point. So these, I start with this as x, and it turns, and it becomes larger into y. In general, in this particular case, all vectors x become longer when they become the y's. But the distinction is two directions, which is the major axis and the minor axis of the ellipse. There, x and y are exactly in the same direction. y is just a different size. So those two directions are called the eigenvectors, or eigendirections of the matrix A. So the eigenvectors of A are two vectors or two directions which don't turn when you multiply them by A. They remain exactly the same. They just get extended or shortened. And the ratio of this length to that length, that ratio is called the eigenvalue. The major axis of the ellipse defines the larger eigenvalue or the dominant eigenvalue. If you take this ratio instead, that's called the second eigenvalue or lambda 2. Now, this is like hundreds of years old problem. So given a matrix, can I find its eigenvalue? That's a lambda 1. So I'm trying to solve. So given A, I'm trying to find lambda 1 and a vector u1 such that A times u1 is lambda 1 times u1. So that's what's happening here. This ratio is lambda 1. And this direction, I don't care about what the magnitude is. The direction is u1. So I can arbitrarily set that the norm of u1 should be equal to 1. It doesn't matter. I started on the unit circle that got stretched into lambda u1. So how do I hunt for the dominant eigenvalue and vector of a matrix? So suppose I initialize some arbitrary x. Two random values, 2d. So x0, x1, initialize to something arbitrary. And then I keep doing the following in a loop. What do I do? I say y equal to ax. And x equals y divided by the norm of y. So anyone remembers what the norm of y is? So norm of y is equal to square root of y0 squared plus y1 squared. We already know how to compute that. We have seen it. So this is the L2 norm. I divide it by the norm. Therefore, the new x again becomes a point on the unit circle. So in effect, going back to the diagram, this is what I'm doing. Let me reuse this. So I start with this x. That maps to the y. I scale down y by going back to the circle. Then I map this to another y. I scale that back and keep going. So I start with this x. I go to the first y. I come back here. I multiply by a. I scale back. Multiply by a. Scale back. So somehow from the picture, it seems like I will be converging on to the dominant eigenvector. But can we prove that that will always happen? So here it goes. So remember I started off with an arbitrary x. Now these two eigenvectors, later on we'll see the major and minor axis of an ellipse are at right angles to each other. So I started off with a coordinate system which looks like this. Those are the two axes. But I can equivalently think of a coordinate system with u1 as the first axis and u2 as the second axis. They're also at right angles. So suppose my initial x could be expressed as a linear combination of the two axes in the new system. So see, I plotted everything in this coordinate system. But I can express any point also as a linear coordinate system in alpha 1 on this axis, alpha 2 on that axis. So any vector I take can also be expressed in the linear combination of u1 and u2. Suppose alpha and alpha 2 are the new coordinates in that system. What happens if I multiply x by a? That is equal to alpha 1 a1 u1 plus alpha 2 a2 u2. But remember by definition, a of u1 is equal to lambda 1 u1. a of u2 equal to lambda 2 u2. So I can replace it here by saying alpha 1 lambda 1 u1 plus alpha 2 lambda 2 u2. This I can write as lambda 1 alpha 1 u1 plus alpha 2 lambda 2 over lambda 1 u2. Now suppose I say a times a times x. That's what happens if we keep on multiplying. This now is equal to lambda 1 alpha 1 a u1 or linear alpha 2 lambda 2 upon lambda 1 a u2. Again by replacing those things I get lambda 1 alpha 1. What is this? This is lambda 1, right? So lambda 1 squared u1 plus alpha 2 lambda 2 over lambda 1 whole squared u2. In general, if I have multiplied x, say k times, I get lambda 1 to the power k alpha 1 u1 plus alpha 2 lambda 2 over lambda 1 to the power k u2. Now remember that lambda 2 is less than lambda 1. Like the minor axis is shorter than the major axis. Therefore, as k increases, this quantity goes to 0. So this whole second thing disappears. And what I get is some power of lambda 1 times the first eigenvalue vector. But remember that in my pseudo code, every step I was dividing by the norm. So this lambda 1 to the power k actually disappears. I'm compensating for it every iteration. And what remains is just the direction u1. So the summary is that any multiple in front of the second eigenvalue fades away geometrically. And only the direction u1 survives. That's why as I iterated this, I kept on turning closer and closer to the major axis. So that's eigenvalues and eigenvectors in 10 minutes. So now let's quickly look at the code. As before, we'll get a little tired of using C++'s very primitive 2D arrays, which can't remember their size, which can't do index checks. And passing C++ arrays into functions is a major headache. Why that is, we'll come to after midterm when you discuss pointers again. But I think you can trust that dealing with native 2D matrices in C++ is too messy. So we'll use this library called Boost. Boost is a very well-known established library which has lots of support for doing different things. The first one we'll use from Boost is the vector and matrix classes. Just like we included iostream, you can install Boost on your system and include Boost's own implementation of vector, which is different from C++'s. In a different namespace, Boost's implementation of a matrix and Boost's implementation of input output for vectors and matrices. So we'll start with that. Now how do you use that? You say using namespace Boost numeric ublas. This is basic linear algebra system or something like that. Now once you use that namespace, you can pick up something called matrix. Just like vector, matrix is a polymorphic container. You can put doubles in it. You can put ints in it. You can put floats in it. So you say I want a matrix of doubles, which is of size 3 by 3. Note the slight difference in syntax. You give brackets here. And now I can fill up m with various things. So I say unsigned i equal to 0 through size 1. This gives me the number of rows. Size 2 gives me the number of columns. And the reason is that there are higher order matrices supported by Boost, which are like 5, 6, any number of dimensions. Now in this particular case, I initialize ij to be 3i plus j. And then I'll print out m, because I'm not in the standard namespace. I'm using Boost's namespace. To be able to use something from STD, I have to declare STD, including STD and L. So this prints out the matrix in one shot. And true to Boost, Boost has libraries for transposing a matrix in one shot. So m will transpose the matrix and print it to standard out in one line. Boost has something called an identity matrix. I don't need to initialize it like we saw. We can just say identity matrix of doubles of size 3 by 3. And I can print that out. I can have a 0 matrix, 3 by 3. I can also have a scalar matrix, where all elements have the single value 20. Scalar matrix internally doesn't take any storage, because you ask for any cell. It will just return 20 to you. And finally, I have a double matrix, where every element is 1.213, double vector. And I take my ms matrix here, and I multiply the matrix divided by 10 and the vector v. So what I'm trying to show is the power of using functions and methods from other people's libraries. I declared a matrix called ms, which has a scalar value, a fixed value of 20 in each cell. I can just say ms by 10, and each cell will become 2. That slash is not c++'s native slash. This slash is a matrix slash, which divides each element of a matrix by a value. And I can now do a matrix vector product by just calling one function from boost, prod of a matrix and a vector. And now I'll print the product of the matrix and the vector. So I'm doing a bunch of things here. So this compiles just like any old function once you include those vital things. Takes longer because it's including the boost library and crunching through all that. Now if I run this, the first thing I did was declare this CYC matrix and initialize with 3 into i plus j. So that is the matrix looked at row y's is 0, 1, 2, 3, 4, 5, 6, 7, 8. That's just a cell number, in other words. Now, transpose of that matrix is what? 0, 3, 6, 1, 4, 7, and 2, 5, 8. That's a transpose. Printed all in one line. No need for coding. Next one is identity matrix of size 3 by 3. 1, 0, 0, 0, 0, 1, 0, 0, 0, 1. Next is the 0 matrix, which is 3 by 3, all 0's. Next I have the constant matrix with every element being 20. So here it is, 20, 20, all the way. Next is this vector v, where all elements are 1.213. And if I take the 20 matrix divided by 10, then each element is 2. And if I multiply that by this fixed thing, it thrice that because it's 3 by 3. So I get a vector which is 7. something. And that's the product being printed. So the quick summary is that I can now declare 2D matrices, initialize them with all kinds of things. This is the initialization code. Notice that the syntax is different. It's not m box i box j. It's bracket open, i comma j bracket close. So very close to mathematical notation. And I can initialize matrices, vectors, multiply them, transpose them. So we'll continue with this next time.