▲ 333 points
back
72 comments
I wish there was a Missing Semester of Linear Algebra course to help people go from "Okay I have a course or two in linear algebra, I know what span, vectors, basis and dimension mean, the formal definition of an inner product space, and I can do Gauss-Jordan elimination, determinants and eigenvalues for small matrices with paper and pencil" to "I have a 100x100 matrix of noisy data from sensors and this research paper I found tells me I can do some fantastic stuff if I compute such-and-such involving eigenvalues or inverses or whatnot. Or maybe I have a process with 1000 states where I know the probability objects move from state i to state j for each pair (i, j) and I want to find the steady state. How do I wrangle numpy into doing what I need?"
MIT has a course called The Missing Semester of Your CS Education [1]. It tells you about practical stuff that you need to know but isn't really taught in classes (shells, version control, build systems, package managers, VM's).
There needs to be something similar for linear algebra, it seems like there's a lot of folk knowledge and a big gap between what typical undergrad courses train you to do and what you encounter in actual practical problems.
(And don't get me started on all the weird linear algebra stuff they have going on in e.g. quantum physics.)
Stanford has ENGR108 [1] based on freely available book "Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares" [2] by Stephen Boyd, with video lectures [3] available. EE263 [4] is sort of a continuation of this at a more advaneced level, it originally also was developed by Boyd and also has video lectures available [5]
[1] https://stanford.edu/class/engr108/
[2] https://web.stanford.edu/~boyd/vmls/
[3] https://www.youtube.com/watch?v=oR6G1MUMveE
[4] https://ee263.stanford.edu/
[5] https://www.youtube.com/playlist?list=PL06960BA52D0DB32B
Sounds like a vanilla numerical math course?
Sounds like EE263 at Stanford. Stephen Boyd’s lectures are on YouTube and so are most of the slides and psets.
check out prof. steve brunton’s youtube channel[0] — it’s light on the numpy side of things, but he’s got a lot of material going from first principles to dealing with huge matrices, compressed sensing, dynamical systems, ML, etc
I don't exactly know what you mean, but at least most points in the blog post were covered in my undergrad numerics courses.
It's a master degree course called Matrix Analysis (e.g Horn and Johnson as text).
Prof Higham has quite the list of papers on the finer points of numerical algorithms going back to his MSc in 1983:
http://www.ma.man.ac.uk/~higham/papers/bibbase.php
including a brief note on the comparing "Top 10 Algorithms in Applied Mathematics" between 2000 and 2016 that may interest some:
https://nhigham.com/2016/03/29/the-top-10-algorithms-in-appl...
It'd be nice to see some meat on the bones and a few ripping yarns about the application end of applied math techniques .. eg: forming an enhanced image from tens (or hundreds) of thousands of multichannel spectral samples using a sensitivity adjusted SVD, and then removing the most common expected background to highlight the anomalies.
It's dry stuff in Linear Algebra, somewhat more exciting when searching for nuclear weapons in a forest or gold in a desert.
He also wrote "Accuracy and Stability of Numerical Algorithms" which is part of the canon of the field.
Why is Newton’s Method not used more often in ML? I know Newton’s method requires the 2nd derivative while Gradient Descent family of algorithms only requires the first, but shouldn’t Newton’s Method be just as straightforward when using autodifferentiation since the process is just a computational graph graph transform using a lookup table?
Re: the first sin, I was told to not invert matrices but never given a satisfying reason. So I wrote a blog post about it:
http://gregorygundersen.com/blog/2020/12/09/matrix-inversion...
The other obvious reason to avoid inverses is that they're only defined for square matrices, whereas LU decomposition works on general rectangular matrices (or rather PLU decomposition,same basic idea).
Invert can also destroy structure you might want to keep around, for example the LU factorization of a banded matrix will still be banded (lapack will do partial pivoting which will increase the bandwidth, but with it'll only double the number of super-diagonals in L), while the inverse is a full matrix.
Most people should pay extra attention to #7.
E.g. Goodfellow et al did even worse than this sin in the Deep Learning book when they claimed the condition number for a square (but not necessarily normal) matrix is defined in terms of eigenvalues. This is false, but nevertheless see 4.2 in https://www.deeplearningbook.org/contents/numerical.html . When I've raised this with people in real life, I typically get some reflexive response that it should be a useful approximation, but as this blog points out, that isn't true either.
Good catch!
Also, please check other posts in Professor Nick Higham's blog. Especially his wonderful and accessible "What is ..." series of articles [0]
Not sure I 100% agree with (2) forming A^TA.
In many real-world use cases A^TA+choleksy is going to be considerably faster than QR on A, and come with few numerical consequences. Even in the numerically challenging cases, pivoted LDL^T on the saddle point system is still going to be faster than the suggested solution of doing QR on A. (Essentially no optimization solver I've ever seen uses QR).
Love the post. I'll take this opportunity to link to a favorite classic linear algebra paper in a similar vein: "Nineteen Dubious Ways to Compute the Exponential of a Matrix" [1]
[1]: https://www.math.purdue.edu/~yipn/543/matrixExp19-I.pdf
Great list thank you. Unfortunately I commit 4 of these on a regular basis.
Why, and which 4?
... same ... :-/
Numerical algorithms have always been a fascination of mine and I spent quite a bit of time studying them. Linear algebra has always seemed to provide some of the most rich content (followed by differential equations, imo). To me linear algebra was so dry when done on paper but suddenly a new world opened up when I could use computers.
I recommend Numerical Linear Algebra by Trefethen and Bau.
Thanks for this excellent link.
Amazing text indeed, this was our textbook in some graduate courses I took in numerical computation. Machine error epsilon still haunts me at night
Trefethen's video lectures are also superb: https://podcasts.ox.ac.uk/series/scientific-computing-dphil-...
The PDFs for Higham's "What is" series are available on Github: https://github.com/higham/what-is
nit: in "5. Not Exploiting Structure in the Matrix", the author says circulant matrices can be solved in O(n log_2 n) operations, where log_2 means "base 2" log. This notation is unnecessary since different bases of logs (as long as the bases are constant and independent of n) only differ in a constant factor, the _2 is insignificant in the big O notation, so it's just O(n log n)
Learned a lot of these lessons when studying statistics as a non-coder trying to implement my own techniques and other peoples research.
Nice list. I would add
8. [edit: oops that's already no 5] Not taking advantage of matrix structure (symmetric, sparse, banded, Toeplitz, ...)
9. Transposing a matrix (Like the inverse A^{-1}, the explicit transpose A^t is often not needed)
IMO 9 is the fault of the language rather than the programmer. Julia solves this one by having types that can represent conjugates, transposes, and adjoints of arbitrary matrices lazily.
> 9. Transposing a matrix
... if you only use the transpose once. If instead it is going to be used multiple times, explicitly computing the transpose can be a huge performance boost.
For dense matrices, it is typically used to exploit memory locality (i.e. to be prefetch- and cache-friendly).
For sparse matrices (your point 8), the advantage can be even more pronounced, sometimes the difference between being able to exploit sparsity, or not.
> 8. Not taking advantage of matrix structure (symmetric, sparse, banded, Toeplitz, ...)
This is the sin no 5 of the article: "5. Not Exploiting Structure in the Matrix"
> 9. Transposing a matrix (Like the inverse A^{-1}, the explicit transpose A^t is often not needed)
Most matrix libraries should make transpose, conjugate, and conjugate transpose just twiddling a bit on its internal representation--BLAS routines should have a parameter on them saying if the input matrix needs to be transposed and/or conjugated before doing an operation.
Materializing the transpose of a matrix is one of the most common and useful operations at the start of a large-scale data processing system.
People who say "transpose is an O(1) operation because it just creates a view" aren't including the important detail of caches and access patterns impact on performance.
Per the first sin, is there any alternative for computing the view matrix from a camera matrix in 3d graphics? This is a case where inversion "feels" appropriate to me.
He is mostly talking about computational linear algebra problems of a large scale type due to large matrices: the "computational intensity" comes from having really large matrices (kxk for k = 100's, 1 000's, 10 000's, 100 000's, 1 000 000's, ...).
In computer graphics, the situation is often different. Usually, you have small matrices (kxk for k=2,3,4); a huge number of vectors; and you want to apply your matrix to all of those vectors. Very often, these matrices have very well known forms and also known well behaved inverses. There isn't really a significant computational cost in computing the inverse (you'll very often write down its formula by hand), and conditioning is usually not an issue (consider a rotation matrix for example or undoing translations with homogeneous coordinates).
You may need an inverse, but you shouldn’t call .inverse() or inverse(A).
Since transformation matrices have simple structure, you can invert them much much faster.
Ex: inverse(R, u) is (R^T, -R^T * u)
> Indeed one would not solve the scalar (n = 1) system 7x = 21 by computing x = 7^-1 × 21 but rather would carry out a division x = 21 / 7
I remember learning in algebra to solve this equation exactly the way he described and said we don't use. You multiply both sides by 1 / 7 which cancels the 7 on the left side, because 1 / 7 is the inverse of 7.
Now I just implicitly divide both sides by 7, but I'm still solving the equation by using the inverse of 7...
But that is sort of the point. In theory, it's the same. But that's not how you should actually compute it.
The cheap way is to compute 21/7.
The hard way is to compute 1/7 (one floating point operation), and then multiply it by 21 (another floating point operation).
With matrices, the discrepancy in work and possibly precision between "solve Ax=b" and "compute x = A^-1 b" can be very large.
I think the difference is that you did the inversion symbolically, not numerically. This is about numeric computation not applying algebraic manipulations.
In linear algebra finding a general inverse of a matrix is harder than solving a specific equation.
ok so what is sqrt(7)x = 21?
well you just 1/sqrt(7) * sqrt(7)x = 1/sqrt(7) * 21
so x = sqrt(7)^2 * 3 / sqrt(7) = sqrt(7) * 3
but you didn't compute the inverse did you, and neither did i. you factored 21 and used the cancellation law (ax = ay => x = y)
I love this topic. Are there any industry fields or sub fields that make heavy use of this content?
Image manipulation, game engines, graphics, finite element analysis, fluid dynamics, signal processing, weather models, etc. It's used everywhere where you have a system of equations to solve numerically.
Most of the industries are dominated by subject matter experts who write code.
Linear algebra is basically the foundation of most serious computation effort. We rank supercomputers essentially by their sustained floating-point computation rate solving large linear equations--a substantial fraction of HPC workloads boil down to "math on large matrices." In the large scale, basically any simulation (whether it be weather modeling or analysis of the stress on a physical object) is solving linear equations. Graphics tends to boil down to doing lots and lots of small matrix multiplications, as does signal or image processing. Optimization problems also tend to boil down to linear algebra as well (via simplex)--so things like operations research, place & route on FPGAs, or even optimal factory ratios in a game like Factorio turn out to need to use linear algebra.
Control theory, optimization
What is a good resource to self-study numerical linear algebra?
There are many good resources! A few, depending on your inclination:
- For intuition: https://www.youtube.com/watch?v=fNk_zzaMoSs
- For rigor: https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010...
- For code: https://codingthematrix.com/
- For numerical/algorithmic details: https://people.maths.ox.ac.uk/trefethen/text.html
fantastic post!