This is a repository for the course 18.06: Linear Algebra at MIT in Spring 2022. See other branches of this repository for previous semesters.
Instructor: Prof. Steven G. Johnson. Course administrator: Yair Shenfeld.
Lectures: MWF11 in 10-250. See recordings and handwritten notes posted online, along with other materials (slides, further reading) posted below.
Exams: 11am in 10-250, on Fridays 2/25, 4/1, and 5/6. Final exam: 9–12 in 10–250 on Wednesday 5/18.
Recitations:
- Gary Choi: T2 in 2-136, T3 in 2-136 (office hours W10 virtual, R11 in-person 2-246C)
- Haushuo Fu: T11 in 2-136, T12 in 2-136 (office hours M4 virtual, W1 in-person 2-238)
- Sergei Korotkikh: T11 in 2-131, T1 in 2-136 (office hours W6 virtual, R6 in-person 2-231D)
- Yair Shenfeld: T9 in 2-132, T10 in 2-132 (office hours T4 virtual, R5 in-person 2-231)
Undergraduate Assistants: Subha Pushpita, Isaac M Lopez, and Gaurav Arya. Email them at 1806sp22_ua ατ mit.edu for 1-on-1 technical help with Julia or other questions that don't work well over Piazza etc. Virtual office hours: M7, T5:30, F1–3.
Resources: Piazza discussion forum, math learning center, TSR^2 study/resource room, pset partners.
This document is a brief summary of what was covered in each 18.06 lecture, along with links and suggestions for further reading. It is not a good substitute for attending lecture, but may provide a useful study guide. (You can also look at the analogous summaries from Fall 2018.)
- course overview/syllabus
- pset 0: due Friday Feb 4 at 11am (submit your solutions on Canvas).
- video: see recordings link above.
Slides giving the syllabus and the "big picture" of what 18.06 is about. Introduction to thinking about matrices as linear operations, not just as "bags of numbers".
Further reading: Strang, chapter 1, and section 8.1 on linear transformations. 3blue1brown has a nice video on matrix multiplication as composition of linear transformations. If you've forgotten the basics of how to multiply matrices by vectors or matrices by matrices, google for some tutorial material online (e.g. Khan academy) and do a quick brush-up.
- handwritten notes: see link above (at beginning)
- Gaussian-elimination Julia notebook
(Started with quick review of matrix–matrix multiplication from the end of lecture 1.)
Gaussian elimination for Ax=b: I started with the grade-school/high-school viewpoint of writing out three equations in three unknowns, adding/subtracting multiples of equations until we were left with one equation in one unknown. Then, I wrote the same equations in matrix form, and renamed this process "Gaussian elimination": we add/subtract multiples of matrix rows to introduce zeros below the diagonal, i.e. to make the matrix upper triangular U. We then do the same row operations to the right hand side b to get a new vector c. Finally, we solve Ux=c for x by working from bottom (1 equation in 1 variable) to top, a process called "backsubstitution".
To do the same operations to both A and b, a useful trick for hand calculations is to augment the matrix with a new column representing the right-hand side, forming [A b] before starting Gaussian elimination.
What comes next? The problem with expressing Gaussian elimination this way, as operations on individual numbers in the matrix, is that it is impossible to follow the process in detail for anything except a very tiny matrix. We need a higher-level "algebraic" way to express the process, both to help us understand it and also to help us to use it (e.g. to perform additional algebraic transformations afterwards). To do this, we want to express the process, not as operations on individual numbers, but as matrix operations.
Rewrote Gaussian elimination in matrix form: we multiply a matrix A on the left by a sequence of lower-triangular "elimination matrices" Eₙ to arrive at an upper-triangular matrix U = EA. To solve Ax=b, we can think of the earlier process as multiplying both sides on the left by E, the linear operator representing the composition (product) of all of the elimination steps, yielding Ux=EAx on the left and c=Eb on the right.
Further reading: Textbook sections 2.1, 2.2, 2.3. Strang lecture 2 video.
Optional Julia Tutorial (Feb 2 @ 5pm): Zoom
A basic overview of the Julia programming environment for numerical computations that we will use in 18.06 for simple computational exploration. This (Zoom-based) tutorial will cover what Julia is and the basics of interaction, scalar/vector/matrix arithmetic, and plotting — we'll be using it as just a "fancy calculator" and no "real programming" will be required.
- Tutorial materials (and links to other resources)
If possible, try to install Julia on your laptop beforehand using the instructions at the above link. Failing that, you can run Julia in the cloud (see instructions above).
- handwritten notes
- Matrix inverse and LU notebook
- pset 0 solutions
- pset 1 (due Feb 11)
Showed that Gaussian elimination can be viewed as LU factorization:
- Gaussian elimination A ⟿ U=EA (without row swaps) can be thought of as A=LU: factorizing A into a product of two simpler (triangular) matrices (L=lower, U=upper). U is the matrix that you normally get when you do elimination by hand, and L (the inverse of the elimination steps L=E⁻¹, a lower-triangular matrix with 1's on the diagonal) is essentially a "record" of the elimination steps.
L is the matrix that "reverses" Gaussian elimination: it tells you how to get A back from L. Despite this, I showed in lecture that L is actually easier to get than E: all you do is make a diagonal matrix of 1's, and then fill in the multipliers from the elimination steps (flipping subtraction to addition) below the diagonal. So, L just requires bookkeeping, and no computation.
Computing U is hard (elimination is a lot of work, even for a computer), but once you have U and L then many things that you might want to do with A become easy.
- For example, suppose you want to solve Ax=b, given A=LU. Write LUx=b=L(Ux), and let y=Ux. Then Ly=b, and we can solve for y by forward-substitution. Given y, we can then solve Ux=y by back-substitution. Both of these steps are easy because the systems are triangular.
- Moreover, solving Ly=b turns out to exactly correspond to applying the elimination steps from A ⟿ U to b. (The 1's on L's diagonal mean that there are no divisions required, either.)
- This means that we can re-use L and U to solve Ax=b for many right-hand sides. In contrast, if you "augment" A with b and then do elimination (A|b)⟶(U|y), you get the same new right-hand side y but you haven't kept a record of the elimination steps, so if you have a new right-hand side you might naively repeat the whole elimination process (hard!) rather than solving Ly=b (easy!).
- More generally, whenever you have A as a product of "simpler" matrices (e.g. triangular, diagonal, …), you can solve Ax=b by a sequence of "simpler" solves.
Introduction to the concept of a matrix inverse more generally as the matrix that reverses the action of a linear operator. Key ideas from the notebook:
- A⁻¹ is the matrix that does the "reverse" of A: A⁻¹(Ax)=x for any x. It also follows that A is the reverse of A⁻¹: A(A⁻¹x)=x for any x, i.e. (A⁻¹)⁻¹=A.
- That is, if Ax=b, then A⁻¹b=x (for any x). (Equivalently, it gives the solution to Ax=b.)
- It only exists for square, nonsingular matrices A. (i.e. an m×m matrix A must give m nonzero pivots when you do elimination.)
- Equivalently, A⁻¹ is the matrix for which A⁻¹A = A⁻¹A = I (the m×m identity matrix).
- I is an identity matrix, the matrix that gives Ix=x for any x or IA=A and AI=A for any A. There are m×m identity matrices for all m, and when we write "I" we usually infer from context how big an I we mean.
In the next lecture (which will start with the end of this notebook), we will look at calculating inverses more generally (although it turns out that this is something that you should almost never do explicitly, even on a computer!).
Further reading: Textbook sections 2.5, 2.6. Strang lecture 4 video and lecture 3 video. See also "The key reason why A = LU" in section 2.6 of the textbook.
- handwritten notes
- Matrix inverse and Gauss–Jordan
Reviewed matrix inverses and key properties thereof.
Went through how to explicitly compute A⁻¹ by solving AA⁻¹ = I. Essentially, this is just solving Ax=b multiple times, where b is each column of I. A common way to organize this for hand calculation (ugh) is the Gauss–Jordan algorithm (on a 3×3 example that can also be found in the textbook): If we do row operations on A to get I, then the same row operations on I give A⁻¹! To carry this out by hand, we augment (A|I), do ordinary Gaussian elimination to get (U|C), and then do elimination "upwards" to get (I|A⁻¹).
Matrix inverses are mainly a conceptual tool that we use to move matrices around symbolically in equations. Once you are through with your algebraic manipulations, you might end up with an expression like A⁻¹b — but when it comes time to actually compute the answer, you should read "A⁻¹b" as "solve Ax=b for x by the best available method".
Further reading: Textbook sections 2.5, 2.6. Strang video lecture 3.
- Handwritten notes
- LU factorization for real
Brief review of previous topics in LU factorization with some more examples in the notebook:
- How the L matrix entries are just the multipliers from Gaussian elimination. No extra work is required!
- How in practice, one rarely "augments" the matrix with the right-hand side. Instead, you compute A=LU, substitute this into Ax=b=LUx, let c=Ux, solve Lc=b, then solve Ux=c. In particular, solving Lc=b is exactly the same as performing the Gaussian-elimination steps on c. (The "augmented" method is a little easier for human bookkeeping, but has essentially no advantage for the computer.)
Some new information about LU to complete the story:
- Given A=LU, you can efficiently solve multiple right-hand sides, or equivalently the matrix equation AX=B.
- How row swaps lead to the factorization PA=LU: in practice, the computer almost always does row swaps, and always gives you a permutation matrix P (or its equivalent).
We apply PA=LU to Ax=b in much the same way as for LU; the only difference is that we have to first apply the permutation P to b.
Permutation matrices P are a great example of a linear operator that is often easier to understand (and more efficient) if you don't write it as a matrix, but instead write it as a "vector" p of the permuted indices 1…n in the new order. Then Px is just x[p] in Julia (and very similarly in Matlab and Numpy): just make a new vector by extracting the components p₁,p₂,… of x.
Further reading: Textbook sections 2.7 (on permutations; we will talk about transposes soon), and 11.1. Strang video lecture 4 and video lecture 5. For 18.06, I don't expect you to know the details of how the permutation P in PA=LU is constructed even though you don't know the permutation in advance … you only need to know how to use PA=LU if it (or something similar) is given to you … but if you are interested this "partial pivoting" algorithm is described in lecture 21 of Numerical Linear Algebra by Trefethen and Bau, or in many other textbooks on numerical linear algebra.
- handwritten notes
- Computational complexity
- Singular matrices and rank
- pset 1 solutions
- pset 2
Complexity of matrix operations: why matrix × vector or backsubstitution scale like n² for n×n matrices, while matrix × matrix or Gaussian elimination (LU factorization) scale like n³. Matrices much bigger than a few thousand square quickly become impractical, and really large problems are only tractable because they have special structure like sparsity.
Began talking about singular matrices. With a relatively simple example, showed that singular matrices (matrices with not enough pivots) A can still have solutions to Ax=b, but only for certain right hand sides. And when they do have solutions, they have infinitely many solutions. As a measure of how singular a matrix is, we talk about its rank r, equal to the number of pivots that we have after elimination. It will turn out that this number r is closely related to the geometry of solutions to singular (and non-square) matrix problems. The solvable right-hand sides live in "r dimensions", a subspace of all possible right-hand sides — in our example, we had r=2 for a 3×3 matrix, so the allowed right-hand sides lived in a 2d plane. And if the solutions exist to an n×n system with rank r, then we will see that the solutions "live" in n–r dimensions, e.g. a 1d line in our example.
To make this precise, we first have to go back and generalize our notion of a "vector" to a more abstract notion called a "vector space", and from this we can define "subspaces", and see how certain key subspaces related to a matrix A tell us the geometry of its solutions.
Further reading: Textbook sections 2.6 ("The cost of elimination") and 11.1. For next time, textbook sections 3.1 and 3.2.
- Handwritten notes
Introduced vector spaces (informally, a set V of anything you can add x±y and multiply by scalars αx) and subspaces (a subset of V such that vector operations stay in the subspace). Examples of vector spaces include real n-component vectors (ℝⁿ, or ℂⁿ for complex numbers), real m×n matrices, functions f(x) (ℝ↦ℝ, real numbers to real numbers). Examples of subspaces includes planes or lines through the origin in ℝ³, or the origin itself. The goal of this is to break vector spaces into smaller pieces that we can still do linear algebra on (hence the need for a subspace, not just any arbitrary subset). Subspaces are especially important to help us understand the solutions (if any) of Ax=b for non-square matrices A.
For an m×n matrix A, introduced two important subspaces.
-
First, the column space C(A): the set {Ax for all x ∈ ℝⁿ}. This is the set of right-hand sides b such that Ax=b is solvable, and is a subspace of the "output space" ℝᵐ (not ℝⁿ unless m=n!). Equivalently, C(A) is all linear combinations of the columns of A, which we call the span of the columns.
-
Second, the null space N(A) (also called the "kernel"): the set {x such that Ax=0} ⊆ ℝⁿ (i.e., a subspace of the "input space" ℝⁿ). Given any solution x to Ax=b, then x+z is also a solution if z ∈ N(A) (i.e. Az=0 ⟹ A(x+z)=Ax+Az=Ax=b).
These are very important subspaces because they tell us a lot about the matrix A and solutions to Ax=b. As a trivial example, if A is an n×n invertible matrix, C(A)=ℝⁿ and N(A)={0}. Conversely, if A is an m×n matrix of zeros, then C(A)={0} and N(A)=ℝⁿ.
Defined a basis for a vector space as a minimal set of vectors (we will later say that they have to be linearly independent) whose span (all linear combinations) produces everything in the space.
We can solve Ux=0 by breaking up x=(p,f) into the pivot variables p and the free variables f, and solve for p from f, and in fact this is an upper-triangular solve that we can do by backsubstitution. This means that we can choose any f we want and p is determined. If we choose the usual basis (1,0,0,…), (0,1,0,…) for f∈ℝⁿ⁻ʳ and solve for p, this gives us the special solutions, a basis for N(A).
Showed that the nullspace is preserved by elimination (row) operations, but that the column space is not. So, to find N(A), we can instead do elimination and find N(U)=N(A) for the upper-triangular form U. Defined the rank as the number of (nonzero) pivots, the pivot columns, and the free columns. We now want to find all possible solutions to Ax=0.
Further reading: Textbook, sections 3.1—3.3; Strang video lecture 6 and lecture 7. Note that Strang's lectures and book emphasize the "reduced row echelon" ("rref") form, which is essentially a bookkeeping trick (similar to Gauss–Jordan for inverses) to do the back-solves for the special solutions all at once. I will not emphasize rref form this semester, but you can use it if you want. (In practical applications, rref form is virtually never used, and for that matter one doesn't actually use elimination at all to find null spaces; instead, one uses something called the SVD that we will cover later.)
- handwritten notes
Went through two examples of finding the special solutions as a basis for the nullspace. They key point is always this: given the free variables, we can easily solve for the corresponding pivot variables.
We can now find the complete solution (i.e., all solutions) to a non-square linear system Ax=b. Elimination turns this into Ux=c. We look for solutions in the form xₚ + xₙ: a particular solution xₚ plus any vector xₙ in N(A) (specified explicitly from the basis). The particular solution xₚ can be any solution to Ax=b, but the simplest one to find is usually to set the free variables to zero. That is, we write the solution xₚ=(p; 0) where p is the unknown values in the pivot rows, setting the other (free) rows to zero, then plug into Uxₚ=c to get p. Since this extracts the part of U that is upper triangular, we can easily solve it. Then we add in anything in N(A).
I used a 3×4 example matrix A (similar in spirit to one in the textbook, section 3.3) that was rank deficient: after elimination, we only had two pivots (rank r=2) in the first two rows, and a whole row of zeros. Furthermore, its pivot columns were the first and third columns, with the second and fourth columns being free — this is possible (albeit unusual)! Showed that we could still get the special solutions by solving for the pivot variables in terms of the free variables = (1,0,…) etcetera, and we still got dimension n-r (= 2) for the null space. When solving Ax=b by elimination into Ux=c, however, we only got a solution x if c was zero in the same rows as U. If the zero third row of U was matched by a zero third row of c, then we got a particular solution as before by setting the free variables to zero. If the zero third row of U was not matched by a zero third row of c, then there is no solution: b was not in the column space C(A). This gives an easy way to check whether the right-hand-side is in the column space. We went through two example right-hand-sides: one with no solutions, and one with infinitely many solutions.
Further reading: Textbook sections 3.3–3.4, lecture 8. As noted in lecture 7, I'm not emphasizing the trick of "rref" form (used extensively by Strang's book and lectures) which makes hand calculation slightly easier, but might push students towards memorizing formulas (which are useless in the long run).
- handwritten notes
- pset 2 solutions
- short pset 3 (due Wed Feb 23)
Bases, dimension, and independence. Earlier, I defined a basis as a minimal set of vectors whose span gives an entire vector space, and the dimension of the space as the size of the basis. Now, we want to think more carefully about the term "minimal". If we have too many vectors in our basis, the problem is that some of the vectors might be redundant (you can get them from the other basis vectors). We now rephrase this as saying that such vectors are linearly dependent: some linear combination (with nonzero multipliers) of them gives the zero vector, and we want every basis to be linearly independent. The dimension of a subspace is still the number of basis vectors.
What does it mean to be linearly independent? Given a set of n vectors {x₁, ⋯, xₙ}, if we matrix a matrix X whose columns are x₁⋯xₙ, then C(X) is precisely the span of x₁⋯xₙ. To check whether the x₁⋯xₙ form a basis for C(X), we need to check whether they are linearly independent. There are three equivalent ways to think about this:
-
We want to make sure that none of x₁⋯xₙ are "redundant": make sure that no xⱼ can be made from a linear combination of the other xᵢ's.
-
Equivalently, we don't want any linear combination of x₁⋯xₙ to give zero unless all the coefficients are zero.
-
Equivalently, we want N(X) = {0}.
In this way, we reduced the concept of independence to thinking about the null space. Since the nullspace is preserved by elimination, it follows that columns of A are dependent/independent if and only if the corresponding columns of R are dependent/independent. By looking at R, we can see by inspection that the pivot columns form a maximal set of independent vectors, and hence are a basis for C(R). Hence, the pivot columns of A (i.e. the columns of A corresponding to the columns of R or U where the pivots appear) are a basis for C(A).
It follows that the dimension of C(A) is exactly rank(A).
Went through four important cases for an m×n matrix A of rank r. (Note that we must have r ≤ m and n: you can't have more pivots than there are rows or columns.)
-
If r=n, then A has full column rank. We must have m ≥ n (it is a "tall" matrix), and N(A)={0} (there are no free columns). Hence, any solution to Ax=b (if it exists at all) must be unique. (If m > n, the problem is "overdetermined": more equations than unknowns.)
-
If r=m, then A has full row rank. We must have n ≥ m (it is a "wide" matrix), and C(A)=ℝᵐ. Ax=b is always solvable (but the solution will not be unique unless m=n). (If m < n the problem is "underdetermined": more unknowns than equations.)
-
If r=m=n, then A is a square invertible matrix. Ax=b is always solvable and the solution x=A⁻¹b is unique.
-
If r < m and r < n, then A is rank deficient. Solutions to Ax=b may not exist and will not be unique if they do exist.
Cases (1)-(3) are called full rank: the rank is as big as possible given the shape of A. In practice, most matrices that one encounters are full rank (this is essentially always true for random matrices). If the matrix is rank deficient, it usually arises from some special structure of the problem (i.e. you usually want to look at where A came from to help you figure out why it is rank deficient, rather than computing the rank etcetera by mindless calculation). (A separate problem is that of matrices that are nearly rank deficient because the pivots are very small, but the right tools to analyze this case won't come up until near the end of the course.)
Further reading: Textbook sections 3.3–3.4, lecture 9.
Reviewed and broadened differential calculus (18.01 and 18.02) from the perspective of 18.06, where we view a derivative f′(x) as a linear operator acting on a small change in the input (dx) to give you the change in the output (df) to first order in dx ("linearized"). This viewpoint makes it easy to generalize derivatives, to scalar-valued functions of vectors where f′(x) is the transposed gradient (∇f)ᵀ, to vector-valued functions of vectors where f′(x) is the Jacobian matrix, and even to matrix-valued functions of matrices like f(x)=A⁻¹ where f′(x) is the linear operator f′(x)[dA]=–A⁻¹dAA⁻¹.
Derivatives viewed as linear approximations have many important applications in science, machine learning, statistics, and engineering. For example, went over the multidimensional Newton algorithm for finding roots f(x)=0 of systems of nonlinear equations. At each step, you just solve a linear system of equations with the Jacobian matrix of f(x), and it converges incredibly rapidly. Gave an example demo where we solved a 2d version of the famous Thomson problem to find the equilibrium position of N repulsive "point charges" constrained to lie on a circle (more generally, a sphere or hypersphere).
Further reading: This material was presented in much greater depth in our 18.S096: Matrix Calculus course in IAP 2022. The viewpoint of derivatives as linear operators (also called Fréchet derivatives) was covered in lectures 1 and 2, Newton's method was covered in lecture 4, and automatic differentiation was covered in lecture 5 — see the posted lecture materials and the further-reading links therein.
- handwritten notes
- pset 3 solutions
- pset 4: due Friday March 4
Briefly discussed the end of lecture 10's slides: another important use of linear approximation via derivatives is for optimization. The gradient gives you the "uphill" direction, so you can maximize/minimize a function by "walking uphill/downhill", leading to a family of algorithms known as gradient ascent/descent, respectively. There are lots of details I won't go into, but the upshot is that you can optimize functions of millions of variables (or more!), which is the key to machine learning (e.g. deep neural nets), large-scale engineering optimization, and more.
Further reading: See part 2 of Matrix Calculus lecture 4 and further-reading thereof.
Reviewed the dot product or inner product of two real column vectors, x⋅y, defined as ∑ᵢxᵢyᵢ. In linear-algebra terms, we write this as x⋅y=xᵀy in terms of the transpose of the vector x: if x is a column vector, xᵀ is a row vector (sometimes more technically called a "dual" vector or "covector"). The length (or norm) of a vector is is the square root of the dot product with itself: ‖x‖=√xᵀx.
The transpose Aᵀ of a linear operator A is defined so that xᵀAy = x⋅(Ay) = (Aᵀx)⋅y = (Aᵀx)ᵀy: transposes move linear operators from one side to the other in an inner product. In consequence, we find for matrices that Aᵀ is constructed by swapping rows with columns in A. Key properties:
- (AB)ᵀ=BᵀAᵀ
- αᵀ = α for scalars. Hence xᵀy = yᵀx (i.e. x⋅y = y⋅x)
- (A⁻¹)ᵀ = (Aᵀ)⁻¹
Orthogonal vectors are those for which xᵀy = 0. Defined the orthogonal complement S⟂ of a subspace S ⊆ V as {x ∈ V such that xᵀy=0 for all y ∈ S}. Combining a basis for S and S⟂ gives a basis for the whole vector space V, so the dimensions of S and S⟂ sum to the dimension of V.
Taking the orthogonal complements of C(A) and N(A) leads us to the four fundamental subspaces for an m×n matrix A of rank r:
- column space C(A) ⊆ ℝᵐ, dimension r
- C(A)⟂ = left nullspace N(Aᵀ) ⊆ ℝᵐ, dimension m-r
- nullspace N(A) ⊆ ℝⁿ, dimension n-r
- N(A)⟂ = row space C(Aᵀ) ⊆ ℝⁿ, dimension r
Further reading: Textbook sections 3.5, 4.1; video lecture 10, video lecture 14. Julia notebook on transposes and orthogonality.
Update: due to the snow emergency, MIT is closed and the exam will be virtual. It will be posted on github at 11am, and you will have until 2pm to submit your solutions on Canvas. (It is the same 1-hour exam that we would have used in person, but you may use all the time you want between 11am and 2pm … just leave yourself time to submit!) The exam is closed-book/closed-notes and must be completed on your own.
Exam 1 will cover the material through lecture 9 and pset 3, including: linear operators, matrix–matrix and matrix–vector operations and interpretations thereof, writing/working with equations in matrix form, solving systems of equations with one or more right-hand sides, Gaussian elimination, back/forward-substitution and triangular matrices, LU factorization and PA=LU, permutation matrices, matrix inverses and Gauss–Jordan, singular matrices, computational costs (which operations are ~ m² or ~ m³ etc and arranging calculations efficiently), rank of a matrix (= number of pivots), vector space & subspaces, null space & special solutions, pivot/free columns column spaces, bases and dimensions of vector spaces, checking whether Ax=b is solvabe, complete solutions to Ax=b (including non-square matrices).
- (Optional) review session: Thursday 2/24 4–5pm in 3-270. Practice problems from review session: spring 2017 problem 2 (solution), fall 2017 problem 1 (solution), spring 2008 problem 4 (solutions), spring 2017 final problem 5(b–c) (solution), fall 2013 problem 3 (solution)
Practice problems: spring 2017 exam 1 problems 1–4 (solutions); fall 2017 exam 1 problems 1-4 (solutions); fall 2017 exam 1 problem 1a,b (solutions); fall 2014 exam 1 problems 1, 2, 3a, 3b, 3c (solutions); fall 2012 exam 1 problems 1–4 (solutions); spring 2008 exam 1 problems 1–4 (solutions); fall 2011 exam 1 problem 1a, 1b, 2 (solutions); spring 2012 exam 1 problems 1–4 (solutions); fall 2013 exam 1 problems 1, 2, 3, 5 (solutions); fall 2009 exam 1 problems 1, 2, 3a–3d, 4 (solutions); spring 2008 exam 1 problems 1, 3, 4 (solutions)
- handwritten notes
Continued discussion of 4 fundamental subspaces.
In particular, since elimination multiplies A on the left, it multiplies Aᵀ on the right by an invertible matrix. Therefore, C(Aᵀ) = C(Uᵀ), and the pivot rows of U (not A!) are a basis for C(Aᵀ). More importantly, this tells us a very non-obvious fact: rank(Aᵀ) = rank(A). (That is, if you did elimination on Aᵀ, you would get the same number of pivots.)
Talked about the viewpoint that solving Ax=b when A has full column rank (i.e. independent columns) is equivalent to a change of basis, or writing b in a "new coordinate system".
However, viewing the columns of A as a basis or "coordinate system", it becomes clear that some bases are nicer than others. If the basis vectors are nearly linearly independent (e.g. nearly parallel), or equivalently if A is nearly singular, then the new coordinate system can be difficult and confusing to use, e.g. tiny vectors can have huge coordinates (coefficients) in the new basis.
A much nicer "coordinate system" is an orthonormal basis: vectors q₁,q₂,…,qₙ that are orthogonal to one another and have length 1. One way of looking at this: to change "coordinates" to an orthonormal basis just involves dot products.
- If you have a non-orthonormal basis a₁,a₂,…, then to write an arbitrary vector b in this basis, i.e. b = a₁x₁ + a₂x₂ + ⋯ with coefficients x₁,x₂,…, you need to solve a linear system Ax=b for x. Hard! (∼m³).
- For an orthonormal basis q₁,q₂,… then to write b = q₁x₁ + q₂x₂ + ⋯ you can just take dot products xᵢ=qᵢᵀb. For example, if you take the dot product q₁ᵀb, then you get x₁ (the coefficient of q₁), because all the other terms have dot product zero.
Next lecture, we will put these q₁,q₂,…,qₙ into a matrix Q and talk about its properties, leading us to the important concept of "orthogonal" or "unitary" matrices.
Further reading: 3blue1brown has a nice video on changes of basis. Textbook sections 3.5, 4.1; video lecture 10, video lecture 14.
- handwritten notes
Orthonormal bases, matrices Q with orthonormal columns q₁,q₂,… (QᵀQ = I):
- Saying that the columns of Q are orthonormal vectors is equivalent to saying QᵀQ = I.
- It follows that ‖Qx‖=‖x‖, and more generally (Qx)ᵀ(Qy) = xᵀy: dot products and lengths are preserved.
- The projection matrix onto C(Q) is just QQᵀ=q₁q₁ᵀ+q₂q₂ᵀ+⋯. (The rank-1 matrix q₁q₁ᵀ is projection onto the line αq₁.) In general, the q component of a vector can be found just by a dot product.
- Similarly, the least-squares solution x̂ minimizing ‖b-Qx‖ is just x̂=Qᵀb.
If a matrix Q with orthonormal columns is square, then it is called orthogonal (or unitary). In this case, QᵀQ=I means that Qᵀ = Q⁻¹. It also follows that QQᵀ = I: a unitary matrix has orthonormal rows and columns.
The way we've been finding a basis for C(A) etcetera is a conceptually nice way to understand the relationships between rank and dimensions. But it isn't a practical tool. Instead, the main tool for this kind of thing, both in practice and in theory, is the singular value decomposition (SVD). The key observation is this: orthonormal bases are ideal, but if you multiply most orthonormal bases by a matrix A the result is no longer orthogonal. Miraculously, it turns out that there is a special orthonormal basis v₁,…,vᵣ of C(Aᵀ) such that Av₁,…,Avᵣ are orthogonal, and this basis leads to the SVD (next time).
To warm up, we saw that Ax for any x can be written as (Av₁v₁ᵀ + ⋯ Avᵣvᵣᵀ)x, where v₁,…,vᵣ is any orthonormal basis of C(Aᵀ). That is, all that matters to A is the "projection" (components) of x along the v₁,…,vᵣ directions, and we can write A=Av₁v₁ᵀ + ⋯ + Avᵣvᵣᵀ as a sum of rank-1 matrices. In general, these Av₁,…,Avᵣ are not orthogonal, however.
As an aside, considered rank-1 matrices abᵀ for a∈ℝᵐ and b∈ℝⁿ: if a and b are nonzero, then abᵀ is an m×n matrix of rank 1. This is easy to see: every column is a multiple of a, so C(abᵀ)=C(b) is 1d, and similarly the row space is the 1d subspace spanned by b. Since the dimension of the column/row space is the rank, then the rank is 1.
Further reading: Strang, section 4.4, and video lecture 17.
- handwritten notes
- SVD introduction
- pset 4 solutions
- pset 5: due March 11
Introduction to the singular value decomposition (SVD) and its application to low rank approximation.
Started talking about orthogonal projection: for any subspace S ⊆ V, any vector b ∈ S can be written as a sum of two vectors, one in S and one in S⟂. The former is the orthogonal projection Pb of b onto S, where P is the projection matrix, and the latter is b–Pb=(I–P)b (the orthogonal projection of b onto S⟂).
Derived the projection matrix P = aaᵀ/aᵀa onto 1d subspaces with a basis vector a. For a normalized vector q=a/‖a‖, this is just qqᵀ. (In the discussion of the SVD, we saw lots of these projections!)
Further reading: Beware that there are several slightly different forms of the SVD; what I've used here is called the "compact" SVD on Wikipedia and is denoted by an "r" subscript in Strang section 7.2; you can expand this to other forms by augmenting U and V with bases for the left and right nullspaces, respectively, and correspondingly adding rows and columns of zeros to Σ. Strang, sections 7.1–7.2, and video lecture 29. Note, however, that the connection of SVD to eigenvalues/eigenvectors of AᵀA is something that we won't cover until later in 18.06. A cool example that uses the SVD to pull out a few key vectors from a big pile of data is the eigenwalker demo. The name for this technique in statistics is principal component analysis (PCA). Strang section 4.2 on orthogonal projection.
Orthogonal projection onto C(A) and other subspaces, and the projection matrix P.
- Projection matrix P = aaᵀ/aᵀa onto 1d subspaces with a basis vector a.
- Projection matrix P = A(AᵀA)⁻¹Aᵀ onto n-dimensional subspaces C(A), where A is m×n with full column rank (rank n).
- Projection onto C(Q), i.e. a subspace with an orthonormal basis is simply QQᵀ. For example, SVD AV=UΣ ⥰ A=A(VVᵀ + (I-VVᵀ))=UΣVᵀ.
- Key properties P²=P, P=Pᵀ, C(P)=C(A), N(P)=C(A)⟂=N(Aᵀ).
- Projection I-P onto the orthogonal complement of C(A), i.e onto N(Aᵀ).
- Equivalence between orthogonal projection and least-squares: minimizing ‖b-Ax‖ is equivalent to minimizing ‖b-y‖ over y∈C(A), and the solution is y=Ax̂=Pb, where AᵀAx̂=Aᵀb.
Further reading: Textbook 4.2; video lecture 15.
- handwritten notes
- Least-square fitting examples
- Global-warming regression example
Introduced the topic of least-square fitting of data to curves. As long as the fitting function is linear in the unknown coefficients x, showed that minimizing the sum of the squares of the errors corresponds to minimizing the norm of the residual ‖b-Ax‖. Went through several examples (see Julia notebooks).
Derived the fact that minimizing ‖b-Ax‖ or ‖b-Ax‖² (least squares) corresponds to orthogonal projection (hence AᵀAx̂=Aᵀb) using either algebra (showing ‖b-Ax‖²≥‖b-Ax̂‖² for any x) or calculus (setting ∇‖b-Ax‖²=0). Also showed that if A is full column rank then the minimum x̂ is unique (i.e. every x≠x̂ gives a larger ‖b-Ax‖).
Further reading: Textbook section 4.3 and video lecture 16. There are many, many books and other materials on linear least-squares fitting, from many different perspectives (e.g. numerical linear algebra, statistics, machine learning…) The brief discussion at the end of the notebook on Runge phenomena and equally spaced vs. Chebyshev points in polynomial fitting was an entry point into approximation theory; if you are interested, the book by Trefethen and accompanying video lectures are a great place to start.
- handwritten notes
- pset 5 solutions
- pset 6: due March 18, 11am
Showed that N(AᵀA)=N(A), and hence AᵀA is invertible whenever A has full column rank, and AᵀAx̂=Aᵀb is always solvable.
Gram-Schmidt orthogonalization: given a non-orthonormal basis a₁,a₂,…, we can convert it to an orthonormal basis that spans the same space. All we do is to take each vector and subtract the projections onto the previous vectors to make them orthogonal, and divide by their lengths to normalize them.
- Conceptually, it is clearest to go directly from a₁,a₂,… to q₁,q₂,…:
- q₁ = a₁ / ‖a₁‖
- q₂ = (a₂ - q₁q₁ᵀa₂) / ‖⋯‖: subtract the projection q₁q₁ᵀa₂ of a₂ onto q₁ to make them orthogonal.
- q₃ = (a₃ - q₁q₁ᵀa₃ - q₂q₂ᵀa₃) / ‖⋯‖: subtract the projections of a₃ onto q₁ and q₂
- etcetera
- For hand calculation, it is perhaps easier to defer square roots to the end: first we compute a basis v₁,v₂,… that is orthogonal but not normalized and then we normalize to q₁=v₁/‖v₁‖, q₂=v₂/‖v₂‖, etcetera at the end.
- v₁ = a₁
- v₂ = a₂ - v₁v₁ᵀa₂/v₁ᵀv₁
- v₃ = a₃ - v₁v₁ᵀa₃/v₁ᵀv₁ - v₂v₂ᵀa₃/v₂ᵀv₂
- etcetera ... note that the vᵢ vectors are still orthogonal, which is why projecting them is still easy, even though they are not normalized to have unit length.
(The process described above and in the book is known as "classical" Gram-Schmidt. In practice, the computer actually uses more sophisticated algorithms. But classical Gram-Schmidt is still a good way to think about the process.
Further reading: Strang, section 4.4, and video lecture 17.
- handwritten notes
Writing Gram–Schmidt in matrix form: it turns out that what we are "really" doing is factoring A=QR, where Q is a matrix with orthonormal columns spanning C(A) and R is an invertible upper-triangular matrix.
Using QR to solve the least-squares problem: given A=QR, the normal equations AᵀAx̂=Aᵀb turn into the triangular system of equations Rx̂=Qᵀb.
Talked about the computational complexity/scaling: for an m×n matrix A, solving least-squares problem by either the normal equations or by QR scale roughly as ~ mn² (dominated by the the cost of forming AᵀA for the normal equations in the former case, or by the cost of Gram–Schmidt or other algorithms for QR factorization in the latter case).
Some practical tips to keep in mind if you ever need to do least-squares or orthogonal basis for real-world problems (not 18.06 exams!) involving data with finite precision:
-
Never explicitly form AᵀA or solve the normal equations: it turns out that this greatly exacerbates the sensitivity to numerical errors (in 18.335, you would learn that it squares the condition number) Instead, use the A=QR factorization and solve Rx̂=Qᵀb. Better yet, just do
A \ b(in Julia or Matlab) or the equivalent in other languages (e.g.numpy.linalg.lstsq), which will use a good algorithm. (Even professionals can get confused about this.) -
Never use Gram–Schmidt for large matrices, which turns out to be notoriously sensitive to small errors if some vectors are nearly parallel. People still compute the "same" QR factorization, just using different methods! There is an improved version called "modified Gram–Schmidt" described in the book, but in practice computers actually use a completely different algorithm called "Householder reflections." You should just use the built-in qr(A)` function in Julia (or similarly other languages), which will do the right thing most of the time.
This is not unusual: there is often a difference between the way we conceptually think about linear algebra and the more sophisticated tricks that are required to make it work well on large matrices of real data.
The QR factorization derived above from Gram–Schmidt is more precisely called the "thin" QR factorization. There is also a related factorization called "full QR", in which Q is augmented by a basis for N(Aᵀ) to make it a square/unitary matrix and R is augmented by appending m–n rows of zeros. See also this blog post by Nick Higham. In practice, most modern numerical libraries (e.g. Julia's qr function or Numpy's linalg.qr) compute the full QR factorization internally, but store Q implicitly as a sequence of reflection operations, so that you can quickly multiply by Q or Qᵀ without having to store the whole Q matrix.
Further reading: Strang, section 4.4, and video lecture 17. Many advanced linear-algebra texts talk about the practical aspects of roundoff errors in QR and least-squares, e.g. Numerical Linear Algebra by Trefethen and Bau (the 18.335 textbook), but this is beyond the scope of 18.06.
Another wonderful and far-reaching application of these ideas is to realize that the same concepts of orthogonal bases and Gram–Schmidt can be applied to any vector space once we define a dot product (giving a so-called Hilbert space, though we won't use that level of abstraction much in 18.06). In particular, it turns out to be especially powerful to think about orthogonal/orthonormal bases of functions. Introduced a dot product f⋅g=∫fg for functions defined on x∈[-1,1], and we'll use it to make an orthogonal basis of polynomials, the Legendre polynomials.
Further reading: Regarding orthogonal polynomials, for example, see these TAMU notes on orthogonal polynomials and these 18.06 notes on orthogonal polynomials, or many other sources, e.g. this book by Szego (1975).
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent