5 ms·

What Is a Schur Decomposition? (2022)

90 points by vector_spaces 3y ago

jjgreen 3y ago
It is a terrible shame that this excellent "What is .." series is concluded, Higham died earlier this year aged 62.
- stabbles 3y ago
  It's chilling to see the blog archive end at "January 2024"
- an_d_rew 3y ago
  Oh no! I had no idea! I met Higham at a conference many years ago, and the lecture he gave was inspiring! From https://ilasic.org/ilas-net-archive-2024/ https://ilasic.org/ilas-net-archive-2024/ 2024 January 22 ILAS-NET Message no. 2515 CONTRIBUTED ANNOUNCEMENT FROM: Francoise Tisseur and Des Higham SUBJECT: Nick Higham (1961--2024) With great sadness we report that Nick Higham, Royal Society Research Professor and Richardson Professor of Applied Mathematics at the University of Manchester, passed away on January 20, 2024, at the age of 62 after an 18 month struggle with a form of blood cancer. An obituary describing Nick's research and leadership contributions will appear in SIAM News in due course. Francoise Tisseur and Des Higham
- deleted 3y ago
  [deleted]
- krosaen 3y ago
  Wow, apparently published his last post, "What Is the Spectral Radius of a Matrix?" just 8 days before his passing. Legend. https://nhigham.com/2024/01/ https://nhigham.com/2024/01/
stabbles 3y ago
For large and sparse matrices you can use the (restarted) Arnoldi method to compute a partial Schur decomposition AQ=QR where Q is tall and skinny and R has a few dominant eigenvalues on the diagonal (i.e. eigenvalues on the boundary of the convex hull). MATLAB uses ARPACK's implementation of this when you call `eigs` I wrote my own implementation ArnoldiMethod.jl in julia, which unlike MATLAB/ARPACK supports arbitrary number types, and also should be more stable in general, and equally fast. [1] https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl
- sfpotter 3y ago
  Why should it be more stable?
  - stabbles 3y ago
    All Arnoldi methods basically repeatedly grow and shrink a search space until they contain good approximations to Schur vectors. Shrinking (=restart) is necessary to keep constant memory requirements, as well as performance. Restart is somewhat non-trivial because it has to preserve the Arnoldi relation. ARPACK uses Implicit Restart, which is a very beautiful trick based on "perfect shifts": a single iteration of the shifted QR algorithm where the shift is taken as an exact eigenvalue is a way to reorder basis vectors (moving vectors you want to remove to the back, so you can truncate them) while preserving the Arnoldi relation. However, perfect shifts in the QR algorithm are good on paper but can lead to arbitrary loss of precision in finite arithmetic, especially with nearly repeated eigenvalues. I think newer implementations (like SLEPc) are often based on Stewart's Krylov-Schur paper. It's not as pretty but certainly more stable. In short: give ARPACK a matrix with a cluster of almost repeated eigenvalues at the boundary of the convex hull, and it tends to fail.
    - sfpotter 3y ago
      Thanks, very helpful. I’ve used ARPACK extensively and have had no complaints with it, but I haven’t tried this particularly pathological case.
- ilayn 3y ago
  I guess it is a rite of passage to rewrite it. I'm doing it for SciPy too together with Propack in [1]. Somebody already mentioned your repo. Thank you for your efforts. [1]: https://github.com/scipy/scipy/issues/18566 https://github.com/scipy/scipy/issues/18566
  - deleted 3y ago
    [deleted]
  - 3y ago
platz 3y ago
explain matrix functions in more detail and how computing a matrix function on an upper triangular matrix is an advantage. how is this a fundamental tenet of numerical linear algebra? also, if I just want the Eigendecomposition, why do I need schur decomposition? The motivation needs more explanation.
- epgui 3y ago
  Well, either it needs more explanation, or the audience is people with certain prerequisite knowledge.
- cat_man 3y ago
  Matrix functions (at least ones I learned about way back when) are Taylor series representations of a function with the matrix plugged in. For example: exp(A) = I + A + A^2/2 + ... + A^n/n! For matrices, you can show these end up as finite series (i.e., it can be written with the highest n being the dimension of the matrix) The A^n factors can be computed with fewer flops than dense matrices, I believe. The "Q's" don't involve little work since they're "unitary" and satisfy QQ = I, so that for example A^2 = (QTQ)^2 = (QTQQT) = QT^2Q* That's where the f(A)= Qf(T)Q* part comes from. It sounds to me like he's quoting what someone described as a "fundamental tenet of numerical linear algebra" as "anything Jordan decomposition can do, Schur can do better". He doesn't seem to be defending "this is a fundamental tenet", but is saying matrix functions illustrate the "Schur can do better" concept. Both Jordan decomposition and Schur are QTQ* decompositions, so I think he's justifying "Schur can do better" with his comments on poorer stability of Jordan form (i.e., the matrix function property depends on A = QTQ*, which both cases satisfy, but Schur has better numerical properties). I don't think the post said you need the Schur decomposition if you want an eigendecomposition (I'm taking to mean the full eigenvalue + eigenvector pair). It just pointed out that the diagonals are eigenvalues, the first column of Q will be have the eigenvector for the first eigenvalue, and you can obtain eigenvectors from the remaining info if you want it. I'm not sure what the typical approach is for doing eigendecompositions, it might not directly involve Schur decompositions at all. So I agree partially that you could add a little motivation to make it more self contained, but also agree it's probably aimed at a less general audience, like students who'd have more context around some stuff where the details are light.
gmfawcett 3y ago
Remember: if you're not certain that your matrix meets the prerequisites for applying the Schur decomposition, you can always apply the Unschur decomposition instead.
- esafak 3y ago
  Found the dad :)
  - marton78 3y ago
    The undad
- jjgreen 3y ago
  Bravo
lxe 3y ago
This needs a "visual explanation" treatment. The concepts are graspable by the average HN reader, but this presentation is just difficult to approach. I tried using ChatGPT to decompose a lot of the terms, which actually helps, but it's hard to vouch for its accuracy: https://chat.openai.com/share/709f61b3-b4cb-48df-b97f-fa6877e4b549 https://chat.openai.com/share/709f61b3-b4cb-48df-b97f-fa6877...

90 points by vector_spaces 3y ago

jjgreen 3y ago

It is a terrible shame that this excellent "What is .." series is concluded, Higham died earlier this year aged 62.

stabbles 3y ago
It's chilling to see the blog archive end at "January 2024"
an_d_rew 3y ago
Oh no! I had no idea! I met Higham at a conference many years ago, and the lecture he gave was inspiring! From https://ilasic.org/ilas-net-archive-2024/ https://ilasic.org/ilas-net-archive-2024/ 2024 January 22 ILAS-NET Message no. 2515 CONTRIBUTED ANNOUNCEMENT FROM: Francoise Tisseur and Des Higham SUBJECT: Nick Higham (1961--2024) With great sadness we report that Nick Higham, Royal Society Research Professor and Richardson Professor of Applied Mathematics at the University of Manchester, passed away on January 20, 2024, at the age of 62 after an 18 month struggle with a form of blood cancer. An obituary describing Nick's research and leadership contributions will appear in SIAM News in due course. Francoise Tisseur and Des Higham
deleted 3y ago
[deleted]
krosaen 3y ago
Wow, apparently published his last post, "What Is the Spectral Radius of a Matrix?" just 8 days before his passing. Legend. https://nhigham.com/2024/01/ https://nhigham.com/2024/01/

stabbles 3y ago

For large and sparse matrices you can use the (restarted) Arnoldi method to compute a partial Schur decomposition AQ=QR where Q is tall and skinny and R has a few dominant eigenvalues on the diagonal (i.e. eigenvalues on the boundary of the convex hull). MATLAB uses ARPACK's implementation of this when you call `eigs` I wrote my own implementation ArnoldiMethod.jl in julia, which unlike MATLAB/ARPACK supports arbitrary number types, and also should be more stable in general, and equally fast. [1] https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl

sfpotter 3y ago
Why should it be more stable?
- stabbles 3y ago
  All Arnoldi methods basically repeatedly grow and shrink a search space until they contain good approximations to Schur vectors. Shrinking (=restart) is necessary to keep constant memory requirements, as well as performance. Restart is somewhat non-trivial because it has to preserve the Arnoldi relation. ARPACK uses Implicit Restart, which is a very beautiful trick based on "perfect shifts": a single iteration of the shifted QR algorithm where the shift is taken as an exact eigenvalue is a way to reorder basis vectors (moving vectors you want to remove to the back, so you can truncate them) while preserving the Arnoldi relation. However, perfect shifts in the QR algorithm are good on paper but can lead to arbitrary loss of precision in finite arithmetic, especially with nearly repeated eigenvalues. I think newer implementations (like SLEPc) are often based on Stewart's Krylov-Schur paper. It's not as pretty but certainly more stable. In short: give ARPACK a matrix with a cluster of almost repeated eigenvalues at the boundary of the convex hull, and it tends to fail.
  - sfpotter 3y ago
    Thanks, very helpful. I’ve used ARPACK extensively and have had no complaints with it, but I haven’t tried this particularly pathological case.
ilayn 3y ago
I guess it is a rite of passage to rewrite it. I'm doing it for SciPy too together with Propack in [1]. Somebody already mentioned your repo. Thank you for your efforts. [1]: https://github.com/scipy/scipy/issues/18566 https://github.com/scipy/scipy/issues/18566
- deleted 3y ago
  [deleted]
- 3y ago

platz 3y ago

explain matrix functions in more detail and how computing a matrix function on an upper triangular matrix is an advantage. how is this a fundamental tenet of numerical linear algebra? also, if I just want the Eigendecomposition, why do I need schur decomposition? The motivation needs more explanation.

epgui 3y ago
Well, either it needs more explanation, or the audience is people with certain prerequisite knowledge.
cat_man 3y ago
Matrix functions (at least ones I learned about way back when) are Taylor series representations of a function with the matrix plugged in. For example: exp(A) = I + A + A^2/2 + ... + A^n/n! For matrices, you can show these end up as finite series (i.e., it can be written with the highest n being the dimension of the matrix) The A^n factors can be computed with fewer flops than dense matrices, I believe. The "Q's" don't involve little work since they're "unitary" and satisfy QQ = I, so that for example A^2 = (QTQ)^2 = (QTQQT) = QT^2Q* That's where the f(A)= Qf(T)Q* part comes from. It sounds to me like he's quoting what someone described as a "fundamental tenet of numerical linear algebra" as "anything Jordan decomposition can do, Schur can do better". He doesn't seem to be defending "this is a fundamental tenet", but is saying matrix functions illustrate the "Schur can do better" concept. Both Jordan decomposition and Schur are QTQ* decompositions, so I think he's justifying "Schur can do better" with his comments on poorer stability of Jordan form (i.e., the matrix function property depends on A = QTQ*, which both cases satisfy, but Schur has better numerical properties). I don't think the post said you need the Schur decomposition if you want an eigendecomposition (I'm taking to mean the full eigenvalue + eigenvector pair). It just pointed out that the diagonals are eigenvalues, the first column of Q will be have the eigenvector for the first eigenvalue, and you can obtain eigenvectors from the remaining info if you want it. I'm not sure what the typical approach is for doing eigendecompositions, it might not directly involve Schur decompositions at all. So I agree partially that you could add a little motivation to make it more self contained, but also agree it's probably aimed at a less general audience, like students who'd have more context around some stuff where the details are light.

gmfawcett 3y ago

Remember: if you're not certain that your matrix meets the prerequisites for applying the Schur decomposition, you can always apply the Unschur decomposition instead.

esafak 3y ago
Found the dad :)
- marton78 3y ago
  The undad
jjgreen 3y ago
Bravo

lxe 3y ago

This needs a "visual explanation" treatment. The concepts are graspable by the average HN reader, but this presentation is just difficult to approach. I tried using ChatGPT to decompose a lot of the terms, which actually helps, but it's hard to vouch for its accuracy: https://chat.openai.com/share/709f61b3-b4cb-48df-b97f-fa6877e4b549 https://chat.openai.com/share/709f61b3-b4cb-48df-b97f-fa6877...