OK, hi guys. Welcome to CS445 Computational Photography. I'm Ashwin Ramesh, and today I'll be doing a brief prerequisite review of some of the linear algebra concepts that you might need to use in this course.
Note that this should just serve as a guide for what you should review before you embark on the projects, rather than a comprehensive review.
With that, let's get started.
So the goals of this session are as follows:
We're going to talk about systems of linear equations and how to represent them as matrices; go over some basic notations, properties,
setting up systems, and concepts like SVD.
And we're also going to do a demo on Jupyter Notebook of some of the library functions in NumPy
that allow you to do these things.
OK. So, here are some references that you should consult.
The YouTube channel 3Blue1Brown in particular is one that I like a lot.
It will really help you get an intuition for linear algebra concepts.
Really great visualizations. I encourage you to check that out.
The first one, the Stanford review PDF, is a good reference for sure.
OK. Let's move on.
So, this is a system of linear equations. This particular system has two equations and two variables.
So, we know it's linear because the highest power of any variable in the equations is 1.
And we know that these...
So, what else do we know? We also know that these equations are linearly independent.
Why do we know this? You can pause and take a look at it, if you would like, but the reason we know is because...
Because the equations are not linear transformations of each other.
So let's look at how to represent this in matrix form. This system is in the form Ax = b,
where A is the matrix of coefficients of the variables and b is a vector of the right hand side constants.
Something to note in this is that you can interchange the equations above,
and that doesn't really affect the solution values of x1 and x2.
So consequently, you can interchange the rows in A and b, and that won't affect the solution to the system.
And, we actually use this property quite extensively when we're solving systems of linear equations using row reduction, and etc.
Another benefit of converting the system into matrix form is that we can employ a lot of parallelism (thread-based parallelism).
Matrices are easy to represent in computers.
OK, so let's look at some basic notations.
So the i-th element of a vector x (and usually the vectors we deal with in this course are column vectors),
is denoted by xi. So the i-th element is denoted by xi,
a matrix "A" with m rows and n columns is said to be in R m cross n, R^(mxn).
And each element is denoted by a lowercase "a" subscripted with the row index, then the column index.
Let's look at some properties of vectors.
A norm of a vector: As you know, a vector is a quantity with direction and magnitude, and a norm is some measure of the magnitude of a vector.
And, a particular norm that's important to remember is the L2 or Euclidean norm. And, it's denoted like that:
It's basically equal to the square root of the sum of the squares of each element in the vector.
So, some vector-vector relationships:
Given two column vectors x and y of the same length, the quantity x transpose y (x^T y) is called the "inner product" or "dot product" of the vectors.
And, it's computed as shown: You basically take corresponding elements, multiply them, and then sum together all the products,
in each vector.
The dot product is a very important operation in machine learning and image processing, etc.,
but more generally, as well.
There's also something called the "outer product", and this operation doesn't really require the input vectors to be of the same length.
So, given an m-length vector and an n-length vector, the result x y transpose (x y^T) is going to be of dimensions m cross n.
And the entries of this matrix, each row is going to be the corresponding value in the first vector multiplied by the second vector,
and each column is consequently the first whole vector multiplied by the corresponding column in the second vector.
So, you can observe that relationship.
OK. So now, let's talk about some operations on matrices.
The main one is matrix multiplication.
The product of two matrices A (with m rows and n columns) and B (with n rows and p columns) is the matrix C (with m rows and p columns).
So the most important thing to note is the number of columns in the left matrix must be equal to the number of rows in the right matrix.
So, matrix multiplication really depends on the order in which you multiply the matrices A and B. (We'll discuss about that in a second.)
And each element of C, Cij, is calculated as shown below: Cij equals the sum from k=1 to n of Aik Bkj.
The way I like to visualize it is that each row of the left matrix A
gets dotted with each column of the right matrix B to produce one element of the matrix C.
That's the best way to look at it, for me.
So let's look at some properties of matrix multiplication. It is associative, which means that if you have three matrices multiplied together,
you can do any combination first between AB and BC there: (AB)C = A(BC).
Matrix multiplication is also distributive: A(B+C) = AB + AC. But it is in general not commutative.
As I mentioned earlier, it can be that AB is not the same as BA.
Sometimes BA is not even legal. It's not even defined. (The dimensions must work out as described earlier.)
Another operation we can do on a matrix is the transpose.
The transpose is basically just an interchange of rows and columns.
The rows become columns, and all the columns become rows.
And the transpose has the following properties:
The transpose of the transpose of A, (A^T)^T, just gives back the original matrix A.
The transpose of the product of two matrices is interesting: (AB)^T = B^T A^T
The left product, AB, the transpose of that is just the right product, B transpose A transpose.
So that's important to note.
And then A plus B, the whole transpose: (A+B)^T, is just A^T + B^T. So the transpose there is distributive.
Matrices have a property known as "rank", and that's given by the number of linearly independent columns and rows of the matrix.
And if the rank of the matrix is equal to the minimum of its number of rows and number of columns then it's said to have full rank.
So the rank of a matrix can be at most the minimum number of rows or columns, and the way you can think about that is,
what's the most linearly independent rows you can have in an m by n matrix?
And what's the most linearly independent columns you can have in an m by n matrix?
Yeah, so I encourage you to think about that.
So now the inverse of a matrix: The inverse is only defined for
square matrices. That's important to note.
And it's denoted by "A inverse", A^(-1), as shown there.
And it's the unique matrix such that A^(-1) A is equal to A A^(-1) and they're both equal to the identity matrix of size n.
So that's important. The identity matrix of size n is an n by n square matrix with all 0s, except 1s on the major diagonal.
OK. Now let's talk about eigenvalues and eigenvectors. Given a square matrix A in R^(m x n),
lambda is an eigenvalue of A with a corresponding eigenvector of x if and only if Ax = lambda x,
and x is not a zero vector, so it's not a trivial solution.
The determinant of A is equal to the product of its eigenvalues,
so the product of all possible values of lambda is the determinant of the matrix.
If you remember, the determinant of the matrix is proportional to this area of
this n-dimensional parallelepiped, composed of the column vectors in A.
It also happens to be the case that the sum of the eigenvalues is equal to the trace of the matrix A, which is the sum of the values on its major diagonal.
So in the case of the example shown below, A has columns [1, 3], [3, 2].
The trace would be 1+2 = 3. And that would be the sum of the eigenvalues.
OK, now we're going to talk about diagonalizable matrices.
So, for a matrix A in R^(n x n), if A has n linearly independent eigenvectors, then A is said to be diagonalizable.
So in math that looks like being able to decompose A into the form U, D, U inverse.
This is a diagonalization of A. We're going to talk more about this, but
basically what the definition says is that a matrix A is diagonalizable if and only if
it has n linearly independent eigenvectors, and here we'll look at why that's the case.
So on the left here we have essentially the definition of
what it means to be an eigenvector, so u1, u2, all the way to un, are the eigenvectors of A.
and lambda: this should be λ1 (lambda 1), λ2, all the way to λn.
These lambda values are the corresponding eigenvalues of A.
Now imagine kind of blocking out or expressing these as their corresponding matrices here. This is just a diagrammatic representation of what that looks like, these equations.
So here you have A, and then what you can do is stack the column vectors,
(stack the eigenvectors as column vectors,) into this matrix, right here,
this n cross n matrix where each column of this matrix is an eigenvector of A.
And on the right-hand side here you have another n by n matrix where each column is going to be λi ui for i from 1 to n.
Essentially what this is showing us is this entire system of linear equations.
It's telling us that Au=λu for all matching values of λ and u, for 1 through n.
OK. So...
You can take this right-hand side matrix, and you can actually decompose it into just a matrix of eigenvectors (here, that's the U matrix again), and a diagonal matrix
that is one where only the diagonals are nonzero,
(only the elements of the major diagonal are nonzero,)
where in each column the position at the major diagonal is occupied by the corresponding eigenvalue.
So this is what that looks like.
Everything is zero except for the values in the major diagonal, and those values are just λ1 through λn.
You can actually try this. If you multiply this out, you'll get another n by n matrix, right,
where the first column is λ1u1, the second column is λ2u2, etc.
all the way to the last column being λn un. That's essentially this matrix here.
So what you have is, A times this matrix of eigenvectors equals that same matrix of eigenvectors times this diagonal matrix.
We're going to call this diagonal matrix D.
And what we can do here is, just right-multiply by the inverse of U to get this expression for A:
A=UDU^(-1). Notice that this is only possible if U is invertible.
That means that U inverse exists, U is a non-singular matrix, which requires that each eigenvector be linearly independent of all the others.
It requires that A has n linearly independent eigenvectors.
So that's how we get that A can be decomposed into
a matrix of eigenvectors times this diagonal matrix of eigenvalues times the inverse of the matrix of eigenvectors
if and only if it has n linearly independent eigenvectors.
OK, lastly we're going to talk about singular value decomposition (SVD).
SVD, singular value decomposition is a way to factorize a matrix into three matrices.
So given a matrix A, we can write A as a product of three matrices U, Σ, and V transpose.
Here U and V are orthogonal matrices and Σ is a diagonal matrix.
Remember that orthogonal matrices are two matrices that have completely orthogonal column vectors
which means that the dot product of any pair of those column vectors is zero.
It's also true for orthogonal matrices that their transpose is also their inverse. You can do a quick example of this to find out why that's the case.
So let's look at how to calculate singular value decompositions.
So what exactly are these U and V matrices?
Let's assume that you can decompose A into something like U Σ V^T. [Note that Σ=Σ^Τ because Σ is a diagonal matrix.]
So A^T would be VΣU^T, and the matrix A^T A would be this:
And we know that because U is orthogonal, its transpose is also its inverse, so these two (U^T U) actually become the identity matrix.
So you end up with (A^T)A is V(Σ^T)Σ(V^T) and we know that Σ is a diagonal matrix so its transpose is itself.
And, we know that therefore Σ^T Σ is going to be just a diagonal matrix of the diagonal values squared.
So we let D=(Σ^T) Σ and there we have (A^T) A = V D (V^T),
which basically looks like a diagonalization of A^T A.
So, what are V and V^T? They're just the eigenvectors of A^T A.
Similarly, if you compute A(A^T), that's going to give you UΣ(V^T)VΣ(U^T).
And V is also an orthogonal matrix, so with similar arguments you can show that U is actually the matrix of eigenvectors of AA^T.
So, here that is. We're doing SVD on a matrix M, which is an (m x n) matrix,
and we have that U's columns are the eigenvectors of MM^T, right.
The diagonals of S, which was sigma (Σ) in our previous notation, are the square root of the eigenvalues of (M^T)M.
So this is because we saw that Σ^T Σ is the D matrix in the previous example, here,
right, and these are called the eigenvalues of M^T M so,
the diagonals of S are just the square root of those values.
So to compute these diagonal values you would just
compute the eigenvalues of M^T M and take the square root.
These σ1, σ2, (...) are called the singular values of M.
The columns of V are the eigenvectors of M^T M.
So note that in the decomposition you actually have V^T, so the rows of V^T are the eigenvectors of M^T M.
Good, so that was singular value decomposition.
Now what we're going to do is look at a quick demo of all the linear algebra concepts we've learned and how to use them in NumPy.
OK. So here we are on Jupyter Notebook, and we're just going to go over what we can use to implement these linear algebra concepts we looked at in Python.
The way we generally want to define vectors is as NumPy arrays.
So here we go, here we have three vectors x, y, and z. They're just NumPy arrays.
We can use np.zeros to define zero vectors. We can just give it a 4 here and it's going to give us a zero vector of size four.
And some operations of vectors like we discussed, like the dot product or inner product,
is just a function called numpy.dot (here, imported as np.dot).
I would consult the documentation for np.dot because it's a really versatile function and it takes care of a wide range of use cases.
I would look at this documentation. So if both a and b are 1-D arrays, it's an inner product.
It can also do matrix multiplication when both the inputs are 2-D arrays and there are some other use cases as well.
So, yeah. So let's take a look at...
Yeah, here we go. So that's how we use np.dot for two vectors.
Next is the outer product, and as we know the outer product takes two vectors x and y, and they can be of different lengths, and it returns a matrix.
If x is of length m and y is of length n it returns an mxn matrix.
So, I'd encourage you to investigate that as well.
Let's look at an example here. Let x equal an np.ones of size five: n = np.ones(5)
And then we'll do, y = np.arange((6))
And then we'll do, np.outer(x,y). OK.
You can pause the video here and kind of think about what this might give us.
As we expected, it's a 5x6 matrix, and each row is
the corresponding element of x times all of y,
and each column is the corresponding element of y times x.
OK. Let's look at the Hadamard product, or element-wise multiplication.
So, between two similarly-shaped arrays (arrays that can be broadcast together),
the star operation is the element-wise multiplication operation
so you can also achieve this with a function called np.multiply
Let's actually recalculate this because our x has changed. They couldn't be broadcast together.
We'll redefine x here (to be of size 4), and then run this again.
Our x*z is [1, 2, 3, 0]. Now let's see what this gives us (with np.multiply). Yep, these two are the same.
I'd also look at the documentation for numpy.multiply. Here it tells us exactly what it is. It multiplies arguments element-wise.
All right, moving on.
NumPy has a module known as linalg that has multiple very useful abstractions for linear algebra.
One specifically is the norm. let's take a look at "numpy.linalg.norm"...
Yeah, so, you pass in the input as the first argument,
and the second argument is the type of norm you'd like to calculate.
So this is the L1 norm, and this would be the L2 norm.
The L2 or Euclidean norm is the one we discussed in class.
So, this would be that.
All right, now we'll move on to matrix multiplication.
So let's define two matrices here, A and B.
A is of shape 2x3 and B is of shape 3x4 so our result should be of shape...
Uh, well, pause the video, figure it out, but...
we do have the answer here, and the result is of shape 2x4.
Notice that we've used "a.dot" here, and so the output of a.dot...
or otherwise, we could put this in as np.dot(a,b)...
Oops, I forgot to run these.
When they're 2D arrays as the input, it does matrix multiplication.
You can also use the "at" symbol (@) for matrix multiplication. That will give you the same thing. [But, np.dot and @ behave differently for vectors.] https://stackoverflow.com/q/54160155
Yeah, this is a simple way to represent matrix multiplication as long as the shapes are matched.
There won't be any assumptions on the part of "np.dot" if you do it this way, so you'll be guaranteed that you're doing what you're trying to do.
The transpose is a property of the matrix. You can just call ".T". Let's take a "b.T".
Right. And we can do a "(a@b).T"... Yep, so...
It (the transpose) will just interchange the rows and columns.
linalg allows you to calculate the matrix rank as well, so let's calculate the rank of a, the rank of a@b, yeah. "np.linalg.matrix_rank(a@b)"
This is an important set of things to note.
Now we're going to get to the higher-level functions in linalg, namely computing the inverse.
"np.linalg.inv" is going to be the inverse.
And what this is computing the inverse of is b b^T (as "b.dot(b.T)")
The SVD can be computed using np.linalg.svd.
If you remember, the SVD of a matrix b decomposes it into a matrix of the eigenvectors of bb^T times a diagonal matrix of its singular values times a matrix of the eigenvectors of b^T b, so...
That's why, just make sure that you know what u, s, and vh are here. That's important to know.
And I would look up the documentation for this, to look at its various arguments.
Another thing we're going to need to do is solve linear systems using NumPy, so take this linear system, for example.
We looked at this at the beginning of the linear algebra slides.
We'll represent this as an array: [[4, -5], [-2, 3]] for the matrix of coefficients, and then [-13, 9] for the constant vector.
And then we'll just use x = np.linalg.solve, and "solve" takes the matrix of coefficients first and then the constant vectors.
And the last thing to note is about broadcasting,
Broadcasting is just a set of rules that you need to follow when you're doing operations on two arrays.
To do an operation between two NumPy arrays they need to be compatible, broadcastable.
So when operating on two arrays, NumPy compares their shapes
element-wise starting from the end and working its way to the beginning.
Two dimensions are compatible when they're equal or one of them is 1.
Let's check out an example here: A is a 4-dimensional array of ones, of shape (3,1,2,4).
B is a 4-dimensional array of ones of shape (1,2,2,1).
So let's just compare these arrays the way NumPy would starting from the end, (from the right,)
and going to the left, and try to figure out whether or not they're broadcastable.
So, are these two equal: 4 and 1? No.
But one of them is 1, so these are compatible.
Are these two equal: 2 and 2? Yes. And here in the remaining two dimensions, we have 1s.
So A and B are broadcastable, so let's multiply A and B and check the shape.
OK, there we go. So, what happened was NumPy broadcasted these arrays together and the result had shape (3,2,2,4).
For each dimension, the result will take the non-unit shape, so for the first from the right it took 4, then 2, 2, 3.
OK. Moving on...
Here we created two arrays, and they're of different shapes. The first one is of shape 2x2, the second one is of shape 2x1,
What will the resulting shape be? We know from our procedure above that it'll be 2x2.
And the way NumPy does the broadcasting is that it will duplicate
this column vector (of B) along the dimension that needs to be matched for this array (A).
So, you'll actually end up with a 2x2 array of 3s and that will get added to A.
So that's very important to note: how exactly NumPy does the broadcasting and what it fills into the newly-formed arrays when it broadcasts.
OK! Well, thank you for watching this linear algebra tutorial.
If you have any questions, as always, please feel free to ask on Piazza.