6 ms·

From 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 de

From 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 ago
What makes it the most fundamental?
- edflsafoiewq 4y ago
  It'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 ago
    But the image and kernel are not the only two properties of a matrix.
    - heinrichhartman 4y ago
      Sure, 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?
heinrichhartman 4y ago
As it turns out, the author does cover the rank-decomposition (under the name "Rank-Revealing Factorization") in his blog as well: https://nhigham.com/2021/05/19/what-is-a-rank-revealing-factorization/ https://nhigham.com/2021/05/19/what-is-a-rank-revealing-fact...
- FabHK 4y ago
  Higham there adds the condition that X, Y be well-conditioned (as they are in the SVD, for example). That's a problem with the Jordan decomposition A = X J X^-1: the X is non-singular in theory, but can be ill-conditioned, and thus the Jordan decomposition is not stable and not really used in practice (as highlighted in the original article).
joppy 4y ago
This decomposition is not used much because it is not unique in any sense, it is not numerically stable, and it is fairly uninteresting: the matrix D can always be taken to consist of a diagonal of 1s of length rank(A) and the rest zeros. It is much more interesting once X and Y are constrained to something like unitary matrices (in which case we get SVD). This “rank decomposition” is a bit interesting algebraically rather than numerically, since then the numerical stability problems disappear. Also, if we take all matrices to be over the integers (and require that “invertible” means the inverse is also an integer matrix), then the rank decomposition is a lot like the Smith normal form of a matrix.
- heinrichhartman 4y ago
  It's not unique, but it can be implemented in a way that is stable. E.g. SVD does this (but is quite overkill for the requirements). It's highly relevant from an algebraic perspective, hence it's curious that it's not covered (at all) in the numeric literature.
  - whatshisface 4y ago
    I could be misunderstanding, but I think you'd get all the algebraic content out of it by computing the dimensions of the image and kernel, then working with them from there on. Why would you want to have a matrix decomposition that separated them? Granted, computing the dimensions of the kernel is not so easy, especially because a pair of vectors can be arbitrarily close without being linearly dependent. No wonder there is no stable way to do it, it technically exceeds the capability of finite-precision numbers. Multiplying a vector of unequal components by most numbers, especially those with non-terminating base-two decimal representations, will produce a vector that is linearly independent when rounded back to finite precision. Clearly then, linear independence on computers has to be considered in the continuous sense in which singular values reveal it, where a very small singular value represents an "almost-kernel," which is the closest thing to a kernel you are likely to find outside of carefully constructed examples or integer matrices.
    - heinrichhartman 4y ago
magnio 4y ago
Isn't that just Gaussian elimination?
FabHK 4y ago
Isn't the Jordan decomposition the most fundamental one for theory?

whatshisface 4y ago

What makes it the most fundamental?

edflsafoiewq 4y ago
It'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 ago
  But the image and kernel are not the only two properties of a matrix.
  - heinrichhartman 4y ago
    Sure, 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?

heinrichhartman 4y ago

As it turns out, the author does cover the rank-decomposition (under the name "Rank-Revealing Factorization") in his blog as well: https://nhigham.com/2021/05/19/what-is-a-rank-revealing-factorization/ https://nhigham.com/2021/05/19/what-is-a-rank-revealing-fact...

FabHK 4y ago
Higham there adds the condition that X, Y be well-conditioned (as they are in the SVD, for example). That's a problem with the Jordan decomposition A = X J X^-1: the X is non-singular in theory, but can be ill-conditioned, and thus the Jordan decomposition is not stable and not really used in practice (as highlighted in the original article).

joppy 4y ago

This decomposition is not used much because it is not unique in any sense, it is not numerically stable, and it is fairly uninteresting: the matrix D can always be taken to consist of a diagonal of 1s of length rank(A) and the rest zeros. It is much more interesting once X and Y are constrained to something like unitary matrices (in which case we get SVD). This “rank decomposition” is a bit interesting algebraically rather than numerically, since then the numerical stability problems disappear. Also, if we take all matrices to be over the integers (and require that “invertible” means the inverse is also an integer matrix), then the rank decomposition is a lot like the Smith normal form of a matrix.

heinrichhartman 4y ago
It's not unique, but it can be implemented in a way that is stable. E.g. SVD does this (but is quite overkill for the requirements). It's highly relevant from an algebraic perspective, hence it's curious that it's not covered (at all) in the numeric literature.
- whatshisface 4y ago
  I could be misunderstanding, but I think you'd get all the algebraic content out of it by computing the dimensions of the image and kernel, then working with them from there on. Why would you want to have a matrix decomposition that separated them? Granted, computing the dimensions of the kernel is not so easy, especially because a pair of vectors can be arbitrarily close without being linearly dependent. No wonder there is no stable way to do it, it technically exceeds the capability of finite-precision numbers. Multiplying a vector of unequal components by most numbers, especially those with non-terminating base-two decimal representations, will produce a vector that is linearly independent when rounded back to finite precision. Clearly then, linear independence on computers has to be considered in the continuous sense in which singular values reveal it, where a very small singular value represents an "almost-kernel," which is the closest thing to a kernel you are likely to find outside of carefully constructed examples or integer matrices.
  - heinrichhartman 4y ago

magnio 4y ago

Isn't that just Gaussian elimination?

FabHK 4y ago

Isn't the Jordan decomposition the most fundamental one for theory?