4 ms·

Introduction to the Conjugate Gradient Method Without Agonizing Pain (1994) [pdf]

71 points by chmike 3y ago

mnw21cam 3y ago
I cited this in my doctoral thesis. I'm not sure the title is accurate though. It doesn't manage to remove the "this bit is magic, just do the exact incantations and it'll work out" feel from it.
gh02t 3y ago
Hah, my professor was using this as a reference in my numerical methods class circa 2010. It's a great one, highly recommended.
trostaft 3y ago
I remember being handed this back when I was taking numerical analysis for the first time. It's an old document, but still useful. IMO the critical pieces of CG that make it a favorable choice for many problems in scientific computing are 1) the fact that it can be performed matrix free 2) its rapid convergence behavior on operators with clusters of eigenvalues (useful for low rank structures) Thet being said, practically speaking, even if I know my operator is positive semi definite, I often find minres out performing cg. There's a nice paper comparing that, "CG versus MINRES: An Empirical Comparison".
- stabbles 3y ago
  "minres outperforming cg" likely depends on the stopping criterion, since different norms are used.
  - trostaft 3y ago
    Yes, to (grossly) summarize the conclusion of the paper: minres can outperform CG on the backwards error ||r_k|| / ||x_k|| whereas CG can outperform minres on the absolute error ||x^* - x_k|| and the energy norm ||x^* - x_k||_A.
fastneutron 3y ago
Shewchuk’s work on mesh generation is nothing short of a masterpiece. I will always direct people to the source of his Triangle code as an example of what good, literate C code should look like. His Berkeley page is here: https://people.eecs.berkeley.edu/~jrs/ https://people.eecs.berkeley.edu/~jrs/
- apengwin 3y ago
  Warmly remember taking his computational geometry in undergrad. It was the most engaging and mathematically creative class I took in college.
hazrmard 3y ago
Thank you! I first came across Conjugate Gradient when reviewing the paper on Neural Ordinary Differential Equations. It was quite challenging to parse through that math. This helps.
aqme28 3y ago
Important to note that this method only works on Hermitian (usually AKA symmetric) and positive-definite matrices, both of which are often pretty big qualifiers.
- dataflow 3y ago
  Did positive definite not imply Hermitian? I feel like I'm forgetting my linear algebra. Edit: Yup, Wikipedia agrees "this condition implies that M is Hermitian"; see their counterexample with a complex vector: https://en.wikipedia.org/wiki/Definite_matrix#Consistency_between_real_and_complex_definitions https://en.wikipedia.org/wiki/Definite_matrix#Consistency_be... Note: Crucially, this is specific to the field of complex numbers (hence the discussion of Hermitian vs. just symmetry). For the field of real numbers, PSD does not imply symmetry, though that's commonly assumed for convenience.
  - necroforest 3y ago
    kind of. you can decompose an arbitrary matrix into symmetric and antisymmetric components: R = S + A. Since A = -A^H (anti-symmetric), for any vector x, <x, Ax> = -<x, Ax> => <x, Ax> = 0. So for any matrix where <x, Rx> > 0, you can add an arbitrary anti-symmetric matrix and keep the same induced quadratic form. So people typically enforce symmetry in their definitions because it is the only part that contributes to the quadratic form and is "nicer" to work with (always diagonalizable, positive eigenvalues, etc.) This should generalize easily to the complex/Hermitian case.
    - dataflow 3y ago
      Thanks! But I think you might've missed a subtlety here: > This should generalize easily to the complex/Hermitian case. This doesn't seem to be true, in that it's actually impossible to have a non-Hermitian matrix C such that x†Cx > 0 over the complex numbers for all x. Whereas over the real numbers, with a matrix R, you can have x'Rx > 0 such that R is asymmetric. The subtlety here is that x itself can be complex in the complex case, which further constraints C to be Hermitian - see the Wikipedia link I posted above. In other words, "complex definiteness" is actually a stronger condition than "real definiteness", even for matrices without an imaginary part.

71 points by chmike 3y ago

I cited this in my doctoral thesis. I'm not sure the title is accurate though. It doesn't manage to remove the "this bit is magic, just do the exact incantations and it'll work out" feel from it.

gh02t 3y ago

Hah, my professor was using this as a reference in my numerical methods class circa 2010. It's a great one, highly recommended.

trostaft 3y ago

I remember being handed this back when I was taking numerical analysis for the first time. It's an old document, but still useful. IMO the critical pieces of CG that make it a favorable choice for many problems in scientific computing are 1) the fact that it can be performed matrix free 2) its rapid convergence behavior on operators with clusters of eigenvalues (useful for low rank structures) Thet being said, practically speaking, even if I know my operator is positive semi definite, I often find minres out performing cg. There's a nice paper comparing that, "CG versus MINRES: An Empirical Comparison".

stabbles 3y ago
"minres outperforming cg" likely depends on the stopping criterion, since different norms are used.
- trostaft 3y ago
  Yes, to (grossly) summarize the conclusion of the paper: minres can outperform CG on the backwards error ||r_k|| / ||x_k|| whereas CG can outperform minres on the absolute error ||x^* - x_k|| and the energy norm ||x^* - x_k||_A.

fastneutron 3y ago

Shewchuk’s work on mesh generation is nothing short of a masterpiece. I will always direct people to the source of his Triangle code as an example of what good, literate C code should look like. His Berkeley page is here: https://people.eecs.berkeley.edu/~jrs/ https://people.eecs.berkeley.edu/~jrs/

apengwin 3y ago
Warmly remember taking his computational geometry in undergrad. It was the most engaging and mathematically creative class I took in college.

hazrmard 3y ago

Thank you! I first came across Conjugate Gradient when reviewing the paper on Neural Ordinary Differential Equations. It was quite challenging to parse through that math. This helps.

aqme28 3y ago

Important to note that this method only works on Hermitian (usually AKA symmetric) and positive-definite matrices, both of which are often pretty big qualifiers.

dataflow 3y ago
Did positive definite not imply Hermitian? I feel like I'm forgetting my linear algebra. Edit: Yup, Wikipedia agrees "this condition implies that M is Hermitian"; see their counterexample with a complex vector: https://en.wikipedia.org/wiki/Definite_matrix#Consistency_between_real_and_complex_definitions https://en.wikipedia.org/wiki/Definite_matrix#Consistency_be... Note: Crucially, this is specific to the field of complex numbers (hence the discussion of Hermitian vs. just symmetry). For the field of real numbers, PSD does not imply symmetry, though that's commonly assumed for convenience.
- necroforest 3y ago
  kind of. you can decompose an arbitrary matrix into symmetric and antisymmetric components: R = S + A. Since A = -A^H (anti-symmetric), for any vector x, <x, Ax> = -<x, Ax> => <x, Ax> = 0. So for any matrix where <x, Rx> > 0, you can add an arbitrary anti-symmetric matrix and keep the same induced quadratic form. So people typically enforce symmetry in their definitions because it is the only part that contributes to the quadratic form and is "nicer" to work with (always diagonalizable, positive eigenvalues, etc.) This should generalize easily to the complex/Hermitian case.
  - dataflow 3y ago
    Thanks! But I think you might've missed a subtlety here: > This should generalize easily to the complex/Hermitian case. This doesn't seem to be true, in that it's actually impossible to have a non-Hermitian matrix C such that x†Cx > 0 over the complex numbers for all x. Whereas over the real numbers, with a matrix R, you can have x'Rx > 0 such that R is asymmetric. The subtlety here is that x itself can be complex in the complex case, which further constraints C to be Hermitian - see the Wikipedia link I posted above. In other words, "complex definiteness" is actually a stronger condition than "real definiteness", even for matrices without an imaginary part.