18 ms·

Seven sins of numerical linear algebra

333 points by jeffreyrogers 4y ago

warinukraine 4y ago
I recommend Numerical Linear Algebra by Trefethen and Bau. Thanks for this excellent link.
- rajman187 4y ago
  Amazing text indeed, this was our textbook in some graduate courses I took in numerical computation. Machine error epsilon still haunts me at night
- hemishom 4y ago
  Trefethen's video lectures are also superb: https://podcasts.ox.ac.uk/series/scientific-computing-dphil-students https://podcasts.ox.ac.uk/series/scientific-computing-dphil-...
dcl 4y ago
Learned a lot of these lessons when studying statistics as a non-coder trying to implement my own techniques and other peoples research.
alchemist1e9 4y ago
Great list thank you. Unfortunately I commit 4 of these on a regular basis.
- fho 4y ago
  ... same ... :-/
- actually_a_dog 4y ago
  Why, and which 4?
  - alchemist1e9 4y ago
    > 1. Inverting a Matrix Possibly only a partial sin since I use Moore-Penrose pseudo inverse and L2 ridge regularization. This sin is being commit as a consequence of committing next sin. > 2. Forming the Cross-Product Matrix A^TA Yes yes this is potential numerically bad. However in practice as long as you are careful with scaling it’s perfectly fine and enable better parallel computations. > 7. Using Eigenvalues to Estimate Conditioning This sin is the only way I actually know so after reading this I realized I need to read more on this. The last of the 4, which obviously is an avoidable sin. > 3. Evaluating Matrix Products in an Inefficient Order I definitely have code that should he changed. It’s also on my todo list now to audit a specific routine that I suspect can be fixed. This one is just stupid it can be avoided.
    - jerf 4y ago
      "7. Using Eigenvalues to Estimate Conditioning" Putting this out to invoke Cunningham's Law... my intuition says that while the article may be right about matrices in the real numbers, using the eigenvalues to check for closeness to singularity may be more valid on the floats, because probably what you're testing for isn't "closeness to singularity" but how close you are to having floating point failures, and that seems at least likely to be heavily correlated. I now sit back and wait for someone to explain why this is wrong while I act like this was an entirely unanticipated result of my post.
      - blt 4y ago
        Why would this be correlated more with eigenvalues than singular values?
onos 4y ago
I love this topic. Are there any industry fields or sub fields that make heavy use of this content?
- garvield 4y ago
  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.
  - singhrac 4y ago
    Significantly used in finance as well, risk estimation is often approximately fitting an MVN to sparse data. I’ve used ideas from this blog many times.
  - commandlinefan 4y ago
    Comes up quite a bit in the AI space as well.
- jcranmer 4y ago
  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.
- blt 4y ago
  Control theory, optimization
LouterDesint 4y ago
What is a good resource to self-study numerical linear algebra?
- gwgundersen 4y ago
  There are many good resources! A few, depending on your inclination: - For intuition: https://www.youtube.com/watch?v=fNk_zzaMoSs https://www.youtube.com/watch?v=fNk_zzaMoSs - For rigor: https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010/ https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010... - For code: https://codingthematrix.com/ https://codingthematrix.com/ - For numerical/algorithmic details: https://people.maths.ox.ac.uk/trefethen/text.html https://people.maths.ox.ac.uk/trefethen/text.html
  - generationP 4y ago
    Rigor and Strang are not two words I'd put in the same sentence. Neil Strickland's notes https://neilstrickland.github.io/linear_maths/ https://neilstrickland.github.io/linear_maths/ (for the matrix point of view) or Jim Hefferon's book https://joshua.smcvt.edu/linearalgebra/index.html https://joshua.smcvt.edu/linearalgebra/index.html (for a vector-spacey treatment) are what comes to my mind when I think "rigor" (along with all sorts of older textbooks like Hoffman/Kunze); Jean Gallier's long set of notes https://www.cis.upenn.edu/~jean/math-deep.pdf https://www.cis.upenn.edu/~jean/math-deep.pdf goes even further in that direction.
- Test0129 4y ago
  The cool thing about numerical linear algebra is, as far as numerical analysis is concerned, it's probably one of the most accessible fields to the would-be citizen scientist. It's been a long time since I've taken linear algebra but if you can find an OCW course, or a really good book I'd start there. You'll want some rigor under your belt. Once you do that, any good numerical analysis text will usually cover the basics (you'll need them) before going into the fascinating world of algorithms as they relate to linear algebra.
- jeffreyrogers 4y ago
  Numerical Linear Algebra by Trefethen and Bau is the resource I've seen recommended most often. I've only read the first few chapters but they are good and useful.
sieste 4y ago
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)
- adgjlsfhk1 4y ago
  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.
  - SuchAnonMuchWow 4y ago
    and in python/numpy, transpose is a O(1) operation: it doesn't actually move any data, just changes the strides used to access data points in the underlying buffer.
  - bee_rider 4y ago
    I think they are talking about explicitly transposing it in memory here. A "transpose operation flag" is pretty a standard feature for linear algebra libraries. If someone is writing a higher level object-oriented linear algebra environment, having the matrices carry around a "transposed" flag that they can pass to BLAS seems like a reasonable thing to do. Since, other than in Julia and Matlab, matrix stuff is usually an add-on, we can't really blame the language for this decision IMO. A downside could be -- usually your underlying tuned library will be BLAS and LAPACK, which don't accept a 'transposed' flag for every single operation. So, from the user point of view, it could be kind of confusing -- "when I transpose a matrix and then go on to multiply, the transpose is free. But if I transpose a matrix and then go on to hit it with QR, it for some reason incurs this weird extra cost -- not where I do the operation, but later, in the QR."
- SuchAnonMuchWow 4y ago
  > 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"
defrost 4y ago
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 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-applied-mathematics/ 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.
- bee_rider 4y ago
  He also wrote "Accuracy and Stability of Numerical Algorithms" which is part of the canon of the field.
- melony 4y ago
  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?
  - phao 4y ago
    With Newton's method, you'll be solving Hx=g (H = hessian of f, g = gradient of f) at each iteration. For large number of variables N, building H is of order N^2 and solving Hx = g is of order N^3 with an usual solver. N^2 and N^3 are really large for an already large N. I believe the reason is as simple as that. It isn't that it is tedious and difficult to write down the formulas or the code to compute H. It's just too costly, computationally speaking. There is also an increased memory cost (possibly having to store H). When people do have ways to go around this problem, they do use Newton's method for large scale problems.
lr1970 4y ago
Also, please check other posts in Professor Nick Higham's blog. Especially his wonderful and accessible "What is ..." series of articles [0] [0] https://nhigham.com/index-of-what-is-articles/ https://nhigham.com/index-of-what-is-articles/
_gabe_ 4y ago
> 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...
- eru 4y ago
  In linear algebra finding a general inverse of a matrix is harder than solving a specific equation.
- flqn 4y ago
  I think the difference is that you did the inversion symbolically, not numerically. This is about numeric computation not applying algebraic manipulations.
- nh23423fefe 4y ago
  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)
- FabHK 4y ago
  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.
gwgundersen 4y ago
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/ http://gregorygundersen.com/blog/2020/12/09/matrix-inversion...
- bee_rider 4y ago
  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.
  - bee_rider 4y ago
    Err, the number of super-diagonals in U obviously, haha, oops.
- Mathnerd314 4y ago
  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).
xdavidliu 4y ago
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)
ThouYS 4y ago
fantastic post!
Test0129 4y ago
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.
nsajko 4y ago
The PDFs for Higham's "What is" series are available on Github: https://github.com/higham/what-is https://github.com/higham/what-is
owlbite 4y ago
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).
rpmuller 4y ago
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 https://www.math.purdue.edu/~yipn/543/matrixExp19-I.pdf
bob1029 4y ago
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.
- phao 4y ago
  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).
  - bob1029 4y ago
    Thank you for the clarification. I didn't realize we were talking about that kind of scale.
- mrjet 4y ago
  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)
Py-o7 4y ago
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 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.
- tomrod 4y ago
  Good catch!
csense 4y ago
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.) [1] https://missing.csail.mit.edu/about/ https://missing.csail.mit.edu/about/
- cschmid 4y ago
  I don't exactly know what you mean, but at least most points in the blog post were covered in my undergrad numerics courses.
- bernulli 4y ago
  Sounds like a vanilla numerical math course?
  - Q6T46nT668w6i3m 4y ago
    I don’t know your age but I expect you’re a bit older! This has changed dramatically in the past 20-30 years and has been a contentious issue for the past few in the computer science education community. Many programs, if they even offer numerical analysis coursework, provide it to supplement regular analysis coursework (rather than computer science coursework).

333 points by jeffreyrogers 4y ago

warinukraine 4y ago

I recommend Numerical Linear Algebra by Trefethen and Bau. Thanks for this excellent link.

rajman187 4y ago
Amazing text indeed, this was our textbook in some graduate courses I took in numerical computation. Machine error epsilon still haunts me at night
hemishom 4y ago
Trefethen's video lectures are also superb: https://podcasts.ox.ac.uk/series/scientific-computing-dphil-students https://podcasts.ox.ac.uk/series/scientific-computing-dphil-...

dcl 4y ago

Learned a lot of these lessons when studying statistics as a non-coder trying to implement my own techniques and other peoples research.

alchemist1e9 4y ago

Great list thank you. Unfortunately I commit 4 of these on a regular basis.

fho 4y ago
... same ... :-/
actually_a_dog 4y ago
Why, and which 4?
- alchemist1e9 4y ago
  > 1. Inverting a Matrix Possibly only a partial sin since I use Moore-Penrose pseudo inverse and L2 ridge regularization. This sin is being commit as a consequence of committing next sin. > 2. Forming the Cross-Product Matrix A^TA Yes yes this is potential numerically bad. However in practice as long as you are careful with scaling it’s perfectly fine and enable better parallel computations. > 7. Using Eigenvalues to Estimate Conditioning This sin is the only way I actually know so after reading this I realized I need to read more on this. The last of the 4, which obviously is an avoidable sin. > 3. Evaluating Matrix Products in an Inefficient Order I definitely have code that should he changed. It’s also on my todo list now to audit a specific routine that I suspect can be fixed. This one is just stupid it can be avoided.
  - jerf 4y ago
    "7. Using Eigenvalues to Estimate Conditioning" Putting this out to invoke Cunningham's Law... my intuition says that while the article may be right about matrices in the real numbers, using the eigenvalues to check for closeness to singularity may be more valid on the floats, because probably what you're testing for isn't "closeness to singularity" but how close you are to having floating point failures, and that seems at least likely to be heavily correlated. I now sit back and wait for someone to explain why this is wrong while I act like this was an entirely unanticipated result of my post.
    - blt 4y ago
      Why would this be correlated more with eigenvalues than singular values?

onos 4y ago

I love this topic. Are there any industry fields or sub fields that make heavy use of this content?

garvield 4y ago
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.
- singhrac 4y ago
  Significantly used in finance as well, risk estimation is often approximately fitting an MVN to sparse data. I’ve used ideas from this blog many times.
- commandlinefan 4y ago
  Comes up quite a bit in the AI space as well.
jcranmer 4y ago
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.
blt 4y ago
Control theory, optimization

LouterDesint 4y ago

What is a good resource to self-study numerical linear algebra?

gwgundersen 4y ago
There are many good resources! A few, depending on your inclination: - For intuition: https://www.youtube.com/watch?v=fNk_zzaMoSs https://www.youtube.com/watch?v=fNk_zzaMoSs - For rigor: https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010/ https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010... - For code: https://codingthematrix.com/ https://codingthematrix.com/ - For numerical/algorithmic details: https://people.maths.ox.ac.uk/trefethen/text.html https://people.maths.ox.ac.uk/trefethen/text.html
- generationP 4y ago
  Rigor and Strang are not two words I'd put in the same sentence. Neil Strickland's notes https://neilstrickland.github.io/linear_maths/ https://neilstrickland.github.io/linear_maths/ (for the matrix point of view) or Jim Hefferon's book https://joshua.smcvt.edu/linearalgebra/index.html https://joshua.smcvt.edu/linearalgebra/index.html (for a vector-spacey treatment) are what comes to my mind when I think "rigor" (along with all sorts of older textbooks like Hoffman/Kunze); Jean Gallier's long set of notes https://www.cis.upenn.edu/~jean/math-deep.pdf https://www.cis.upenn.edu/~jean/math-deep.pdf goes even further in that direction.
Test0129 4y ago
The cool thing about numerical linear algebra is, as far as numerical analysis is concerned, it's probably one of the most accessible fields to the would-be citizen scientist. It's been a long time since I've taken linear algebra but if you can find an OCW course, or a really good book I'd start there. You'll want some rigor under your belt. Once you do that, any good numerical analysis text will usually cover the basics (you'll need them) before going into the fascinating world of algorithms as they relate to linear algebra.
jeffreyrogers 4y ago
Numerical Linear Algebra by Trefethen and Bau is the resource I've seen recommended most often. I've only read the first few chapters but they are good and useful.

sieste 4y ago

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)

adgjlsfhk1 4y ago
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.
- SuchAnonMuchWow 4y ago
  and in python/numpy, transpose is a O(1) operation: it doesn't actually move any data, just changes the strides used to access data points in the underlying buffer.
- bee_rider 4y ago
  I think they are talking about explicitly transposing it in memory here. A "transpose operation flag" is pretty a standard feature for linear algebra libraries. If someone is writing a higher level object-oriented linear algebra environment, having the matrices carry around a "transposed" flag that they can pass to BLAS seems like a reasonable thing to do. Since, other than in Julia and Matlab, matrix stuff is usually an add-on, we can't really blame the language for this decision IMO. A downside could be -- usually your underlying tuned library will be BLAS and LAPACK, which don't accept a 'transposed' flag for every single operation. So, from the user point of view, it could be kind of confusing -- "when I transpose a matrix and then go on to multiply, the transpose is free. But if I transpose a matrix and then go on to hit it with QR, it for some reason incurs this weird extra cost -- not where I do the operation, but later, in the QR."
SuchAnonMuchWow 4y ago
> 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"

defrost 4y ago

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 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-applied-mathematics/ 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.

bee_rider 4y ago
He also wrote "Accuracy and Stability of Numerical Algorithms" which is part of the canon of the field.
melony 4y ago
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?
- phao 4y ago
  With Newton's method, you'll be solving Hx=g (H = hessian of f, g = gradient of f) at each iteration. For large number of variables N, building H is of order N^2 and solving Hx = g is of order N^3 with an usual solver. N^2 and N^3 are really large for an already large N. I believe the reason is as simple as that. It isn't that it is tedious and difficult to write down the formulas or the code to compute H. It's just too costly, computationally speaking. There is also an increased memory cost (possibly having to store H). When people do have ways to go around this problem, they do use Newton's method for large scale problems.

lr1970 4y ago

Also, please check other posts in Professor Nick Higham's blog. Especially his wonderful and accessible "What is ..." series of articles [0] [0] https://nhigham.com/index-of-what-is-articles/ https://nhigham.com/index-of-what-is-articles/

_gabe_ 4y ago

> 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...

eru 4y ago
In linear algebra finding a general inverse of a matrix is harder than solving a specific equation.
flqn 4y ago
I think the difference is that you did the inversion symbolically, not numerically. This is about numeric computation not applying algebraic manipulations.
nh23423fefe 4y ago
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)
FabHK 4y ago
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.

gwgundersen 4y ago

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/ http://gregorygundersen.com/blog/2020/12/09/matrix-inversion...

bee_rider 4y ago
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.
- bee_rider 4y ago
  Err, the number of super-diagonals in U obviously, haha, oops.
Mathnerd314 4y ago
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).

xdavidliu 4y ago

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)

ThouYS 4y ago

fantastic post!

Test0129 4y ago

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.

nsajko 4y ago

The PDFs for Higham's "What is" series are available on Github: https://github.com/higham/what-is https://github.com/higham/what-is

owlbite 4y ago

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).

rpmuller 4y ago

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 https://www.math.purdue.edu/~yipn/543/matrixExp19-I.pdf

bob1029 4y ago

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.

phao 4y ago
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).
- bob1029 4y ago
  Thank you for the clarification. I didn't realize we were talking about that kind of scale.
mrjet 4y ago
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)

Py-o7 4y ago

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 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.

tomrod 4y ago
Good catch!

csense 4y ago

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.) [1] https://missing.csail.mit.edu/about/ https://missing.csail.mit.edu/about/

cschmid 4y ago
I don't exactly know what you mean, but at least most points in the blog post were covered in my undergrad numerics courses.
bernulli 4y ago
Sounds like a vanilla numerical math course?
- Q6T46nT668w6i3m 4y ago
  I don’t know your age but I expect you’re a bit older! This has changed dramatically in the past 20-30 years and has been a contentious issue for the past few in the computer science education community. Many programs, if they even offer numerical analysis coursework, provide it to supplement regular analysis coursework (rather than computer science coursework).