9 ms·
The big six matrix factorizations
- isaacfrond 4y agoThe article presents a note on the 6 well known matrix compositions. He states that all of them have cubic complexity, but practical algorithms with better exponents exist for all of them.
- adgjlsfhk1 4y agoThey're all basically O(M(n)) where M(n) is your matrix multiplication time. Even though M(n)<=n^2.3...., it's reasonable to say that it's n^3, because in practice, no one uses the sub-cubic algorithms. Strassen is possibly workable, but it isn't widely used, and all of the sub-cubic algorithms have accuracy tradeoffs.
- photon-torpedo 4y ago> but practical algorithms with better exponents exist for all of them. I'm aware of randomized algorithms with better complexity, which come at the cost of only giving approximate results (though the approximation may be perfectly good for practical purposes). See e.g. [1]. Are there other approaches? [1] https://doi.org/10.1137/090771806 https://doi.org/10.1137/090771806
- owlbite 4y agoIf the data is Sparse (which is not uncommon for large matrices in the real world), you can exploit the sparsity to do significantly better then O(n**3).
- oh_my_goodness 4y agoCould you link a practical algorithm with an exponent lower than 3? (I think of these things https://en.wikipedia.org/wiki/Computational_complexity_of_matrix_multiplication https://en.wikipedia.org/wiki/Computational_complexity_of_ma... as not being practical, but I'd love to be wrong. )
- daniel-cussen 4y agohttps://www.fgemm.com https://www.fgemm.com, coming soon.
- gnufx 4y agoFor something Strassen-ish you could look at https://jianyuhuang.com/papers/sc16.pdf https://jianyuhuang.com/papers/sc16.pdf and the GPU implementation https://apps.cs.utexas.edu/apps/sites/default/files/tech_reports/GPUStrassen.pdf https://apps.cs.utexas.edu/apps/sites/default/files/tech_rep...
- martinhath 4y agoThe fact that they're all cubic isn't really the notable part of the runtime of computing the different decompositions, because the constants involved are really different. In practice, a common reason for computing many of these decompositions is to solve a linear system `Ax=b`, because with the decomposition in hand it is really easy to solve the whole system (using e.g. backsubstitution). For instance, with C++s Eigen, look at the 100x100 column of [1], and we can see that there's orders of magnitude difference between the fast and slow approaches. THey're all still cubic, sure, but we're talking 168x difference here. (of course, it's not so clear cut, since robustness varies, not all methods are applicable, and the benchmark is for solving the system, not computing the decomposition, but overall, knowledge of which decomposition is fast and which is not is absolutely crucial to practitioners) [1]: https://eigen.tuxfamily.org/dox/group__DenseDecompositionBenchmark.html https://eigen.tuxfamily.org/dox/group__DenseDecompositionBen...
- FabHK 4y agoCan you provide a source? For example, the SVD (Golub-Kahan-Reinsch) method generally involves bidiagonalisation of the matrix and then implicit QR steps, which are often achieved with Householder reflections and Givens rotations. Sure, all of those are conceptually matrix multiplications, but they're not implemented as O(N^3) matrix multiplications; rather, their special structure is exploited so that they're faster. Yet, the entire algorithm is still cubic. So not sure Strassen would accelerate things.
- leephillips 4y agoThank you for this wonderfully concise summary: it’s convenient to have all this in one compact document. I suppose “flops” means “floating-point operations” here? Heretofore I’ve always encountered this as an abbreviation for “floating-point operations per second”.
- hanche 4y agoIndeed, the meaning of “flops” is ambiguous, but that seems hard to avoid. Fortunately, the ambiguity is easily resolved from context in most cases, as it is here.
- the_svd_doctor 4y agoIt's always been very ambiguous in practice. To avoid ambiguity I usually use `flop/s` but not everyone likes that :)
- arinlen 4y ago> To avoid ambiguity I usually use `flop/s` but not everyone likes that :) Flop/s makes no sense and is outright wrong. The whole point of flops is to express how many floating point operations are required by an algorithm. How many operations are performed per second is a property of the hardware you're using to run an implementation of the algorithm. You want to express computational complexity in terms of floating point operations.
- zodiac 4y agoMaybe the_svd_doctor uses flop/s to describe hardware, like you said
- gnufx 4y agoThe number of operations per second varies by at least an order of magnitude on the same hardware depending on the GEMM algorithm you use (reference or Goto), and you quote the performance of the implementation in terms of FLOP/s knowing the number of FLOPs required by the computation. That makes sense to people who implement and measure these things.
- melling 4y agoAny suggestions on what to learn in Linear Algebra after Gilbert Strang’s 18.06SC? https://ocw.mit.edu/courses/18-06sc-linear-algebra-fall-2011/ https://ocw.mit.edu/courses/18-06sc-linear-algebra-fall-2011... My goal is to learn the math behind machine learning.
- deleted 4y ago[deleted]
- isaacfrond 4y agoHere are some links from my favourite reference website: Book: Mathematics for Machine Learning (mml-book.github.io) https://news.ycombinator.com/item?id=16750789 https://news.ycombinator.com/item?id=16750789 Mathematics for Machine Learning [pdf] (mml-book.com) https://news.ycombinator.com/item?id=21293132 https://news.ycombinator.com/item?id=21293132 Mathematics of Machine Learning (2016) (ibm.com) https://news.ycombinator.com/item?id=15146746 https://news.ycombinator.com/item?id=15146746
- natly 4y agoThese lectures are fantastic after you've mastered the basics of linear algebra https://www.youtube.com/watch?v=McLq1hEq3UY https://www.youtube.com/watch?v=McLq1hEq3UY (convex optimization, by a very experienced and often funny lecturer)
- nightski 4y agoStephen Boyd is a very good lecturer! I watched his videos on linear dynamical systems almost a decade ago and thought he did a fantastic job. Would highly recommend.
- sjburt 4y agoThe problem set for the stanford linear dynamic systems that he taught is also fantastic. The level of difficulty and the breadth of applications leaves you with so many tools.
- 4y ago
- shusaku 4y agoI love Dr. Higham’s blog. In case anyone missed it, the whole series is on GitHub: https://github.com/higham/what-is https://github.com/higham/what-is
- heinrichhartman 4y agoFrom a theoretical perspective the most fundamental decomposition is the rank-decomposition: A = X D Y with X,Y invertible and D diagonal. It's called rank decomposition, because you can read off the rank from A by counting the non-zero entries in D. It's also useful to determining bases for the image and the kernel of A. Every math student learns a version of that in their first lecture series of Linear Algebra. Curiously, the rank decomposition is not covered in the numerical literature. Also it's not possible to derive a rank decomposition from LR, QR decomposition, although the underlying algorithms (Gauss, Gram-Schmidt) could be used to do so. It took me multiple weeks of work, to understand what the problem with this algorithms is, and what the practical options are to establish a rank decomposition are. Full details are available on my blog: https://www.heinrichhartmann.com/posts/2021-03-08-rank-decomposition/ https://www.heinrichhartmann.com/posts/2021-03-08-rank-decom... TL;DR. Your options are (1) SVD, (2) QR factorization with column pivoting (part of LAPACK), (3) LDU factorization with total pivoting (not implemented in BLAS/LAPACK/etc.)
- whatshisface 4y agoWhat makes it the most fundamental?
- edflsafoiewq 4y agoIt's basically the matrix form of dim(domain(f)) = dim(ker(f)) + dim(im(f)), or domain(f)/ker(f) ≅ im(f) with f (inducing) the isomorphism.
- whatshisface 4y agoBut the image and kernel are not the only two properties of a matrix.
- heinrichhartman 4y agoSure, but from an algebraic (categorical) perspective, taking kernels, co-kernels and images are the most fundamental operations you can think of. The question is: How do you compute them in an effective (avoiding full SVD) and stable way?
- swframe2 4y agoSee this playlist on semidefinite programming: https://www.youtube.com/playlist?list=PLqwozWPBo-FtOsjrla4jk6GYcxX1aN_1a https://www.youtube.com/playlist?list=PLqwozWPBo-FtOsjrla4jk... for the surprising way factorization can be used to partially solve a NP problem.
- howlin 4y agoNow apply them to the data used to identify the "Big Five" personality traits! This is an interesting application of factor analysis/matrix factorization: https://en.wikipedia.org/wiki/Big_Five_personality_traits https://en.wikipedia.org/wiki/Big_Five_personality_traits
- jetbooster 4y agoAnd then sell it to the Big Four! https://en.m.wikipedia.org/wiki/Big_Four_accounting_firms https://en.m.wikipedia.org/wiki/Big_Four_accounting_firms
- Royi 4y agoWhere's the data? It could be indeed interesting.
- FabHK 4y agoYes, the Big Five were found by factor analysis, generally the PCA, which is basically an eigenvalue (spectral) decomposition of the (symmetric positive definite) covariance matrix of the data, which coincides with the SVD of the (normalised) data matrix. The SVD is everywhere. IIRC, the first Netflix price was one won with the SVD: https://en.wikipedia.org/wiki/Netflix_Prize https://en.wikipedia.org/wiki/Netflix_Prize https://cseweb.ucsd.edu/classes/fa17/cse291-b/reading/ProgressPrize2008_BellKor.pdf https://cseweb.ucsd.edu/classes/fa17/cse291-b/reading/Progre...
- whatever1 4y agoBoo LU!!! Go with stable QR!
- the_svd_doctor 4y agoBut LU, when it works, is less flops (only half though, IIRC). And doesn't require a square-root. I think that's why it used to be preferred when compute was a very scarce resource. Clearly that's not really the case anymore today. But otherwise I agree, QR is almost better on every front.
- jdeaton 4y agoThis is a fantastically clear outline of this topic. Thank you! > The terms “factorization” and “decomposition” are synonymous and it is a matter of convention which is used. Our list comprises three factorization and three decompositions. I can't tell if this is a joke: right after saying that these two words mean the same thing in this context, they are then used to categorize the methods. Edit: This is the kind of content that proves the "dead internet theory" is wrong.
- tomnipotent 4y ago> I can't tell if this is a joke It makes sense, the author is just using the original name of the method to bucket them e.g. spectral decomposition.
- yuppiemephisto 4y agothe words are used synonymously in practice, but i think there's an opportunity to distinguish them usefully. suggestion: - 'factoring' is a multiplicative breakdown, a 'composition', like prime factorization. - 'decomposition' could also be called 'partition', and is an additive breakdown, like how 3 could be split into 2 + 1
- FabHK 4y agoIt's not a categorisation, but some trivia regarding the name. For example: By convention, the SVD is called the SVD (singular value decomposition), not the SVF (singular value factorisation). That would be synonymous, but that's not what it's called. Similarly for the other 5 methods.
- akomtu 4y agoImo, a better way to present this is to draw a diagram where various decompositions are connected by arrows that explain how one decomposition can be turned into another.
- jdeaton 4y agoThat is fascinating. Do you have an example of this which you can point to?
- rhyn00 4y agoSpectral decomposition is pretty cool. Application 1 - spectral clustering - an alternative to k-means for nonlinear clusters. Get a Distance matrix of your data, spectral decomp, run k-means on your k top eigen vectors and that's your clusters. Application 2 - graph clustering - (run spectral clustering on adj matrix!) There's some tricks to getting it to work in practice like normalizing but it's a simple and powerful method. Also the matrices can get big so it helps a lot to use sparse matrix libraries for the computations. [1] https://towardsdatascience.com/spectral-clustering-aba2640c0d5b#:~:text=Spectral%20clustering%20is%20a%20technique,non%20graph%20data%20as%20well https://towardsdatascience.com/spectral-clustering-aba2640c0.... [2] https://www.hindawi.com/journals/ddns/2020/4540302/ https://www.hindawi.com/journals/ddns/2020/4540302/
- bloqs 4y agoI love that on this website I can find comments that I simply don't understand a word of. Good on you for doing the stuff you do sir.
- arbitrandomuser 4y agoI've never given a second thought about what the etymology of "spectral" in spectral decomposition is. Somewhere in the back of my mind (and I guess many students of physics have the same notion) subconsciously i assumed it originates from eigenvalues of the Hamiltonian determining the atomic spectral lines . But I've never followed up on it and actually looked it up .
- the_svd_doctor 4y agoI might be wrong about the exact historical reason. But the way I see it "spectral decomposition of A" is a way to express A as a sum of orthogonal, rank-1, operators. A = \sum l_i u_i u_i^T. Those l_i are the eigenvalues; u_i are the eigenvectors. The eigenvectors look a whole lot like the "modes" in a Fourier decomposition. And if you plot (i, l_i), the eigenvalues are a bit like the "spectrum" (the amplitude of each mode). In fact, the complex exponentials (the modes in the Fourier decomposition) are also eigenvectors of a specific operator (the Laplacian). Math people are good at finding connections between things.
- mgaunard 4y agoMissing the LDLT decomposition.
- the_svd_doctor 4y agoIt can't always be computed accurately; I suspect that's why it's not mentioned. Bunch–Kaufman is the "right" factorization for indefinite Hermitian/symmetric matrices but it's not as well know.
- the_svd_doctor 4y agoHigham is an outstanding researcher. One book I really enjoyed from him is https://epubs.siam.org/doi/book/10.1137/1.9780898717778 https://epubs.siam.org/doi/book/10.1137/1.9780898717778. He has many many other excellent ones.
- deleted 4y ago[deleted]
- ted_dunning 4y agoIt is important to mention the role of randomized algorithms for approximate solutions. See this excellent and practical survey. https://arxiv.org/abs/2002.01387 https://arxiv.org/abs/2002.01387
- signaru 4y agoQR factorization leads to many useful stuff : Eigenvalues/vectors, SVD (singular value decomposition), PCA (principal component analysis)
- the_svd_doctor 4y agoQR iteration (which uses QR decomposition) truly is a remarkable and beautiful algorithm. It was even coined one of the top 10 algorithm of the 20th century https://archive.siam.org/pdf/news/637.pdf https://archive.siam.org/pdf/news/637.pdf