4
$\begingroup$

I have to write a little finite elements code in C.

I was asked to implement the conjugate gradients method, which I have done. Now, I am looking to improve further the efficiency of my program by using the Incomplete Cholesky Decomposition.

As I understand it, instead of solving the $ Ax = B $ system, I have to solve the $ M^{-1} Ax = M^{-1}B$ system where matrix $M$ is the incomplete Cholesky decomposition of matrix $A$.

However, when I tried to implement this, I found that computing $M$ and then inverting it was more expensive in time than simply using the CG method. I suppose I must have done something wrong, but I really don't see how PCG with Cholesky is implemented.

Improve this question
$\endgroup$
8
  • 4
    $\begingroup$ Note that you don't ever want to compute $M^{-1}$ explicitly. Rather, turn each expression involving $M^{-1}$ times a vector into the solution of a system of equations with $M$ as the matrix. $\endgroup$
    – Brian Borchers
    Commented May 12, 2018 at 20:30
  • $\begingroup$ Also, that's not how you implement a preconditioned CG method. Have you even looked at the wikipedia page? (To reiterate Brian's comment: You'd solve $Mz_{k+1}=r_{k+1}$ using backward and forward substitution based on the incomplete Cholesky factors, not invert and apply anything.) $\endgroup$
    – Christian Clason
    Commented May 12, 2018 at 20:54
  • $\begingroup$ I am not sure to understand what you mean/how to do it in practice. Moreover, computing M (using Cholesky's incomplete decomposition) takes as much time as directly solving $Ax = B$ when I tried it. $\endgroup$
    – Gilles Poncelet
    Commented May 12, 2018 at 20:55
  • $\begingroup$ I have looked up the Wikipedia page yes. However, my main problem is computing M and finding some way to invert it in a cheap and efficient way. $\endgroup$
    – Gilles Poncelet
    Commented May 12, 2018 at 20:57
  • 2
    $\begingroup$ You never compute $M$ directly, that's the whole point -- you compute an incomplete Cholesky factorization of $A$, which you can invert (factor by factor) very cheaply. That's your $M^{-1}$, if you will. Of course, the threshold value for the fill-in is a parameter you have to tune for your specific problem -- if you choose poorly, you indeed won't see any speed-up. $\endgroup$
    – Christian Clason
    Commented May 12, 2018 at 20:58

1 Answer 1

Reset to default
8
$\begingroup$

While preconditioned CG with incomplete Cholesky (ICC) is reasonably straightforward to formulate mathematically, writing an efficient implementation is a non-trivial matter.

Here's some of the things that would make for an efficient implementation:

  • the ICC implementation exploits the sparsity of your matrix $A$ (you'll want to store $A$ and $M$ in some sparse matrix format)
  • $M$ is computationally "easy" to invert (you won't compute the actual inverse $M^{-1}$, just its action on a vector)
  • you set appropriate hyper-parameters for the ICC algorithm (e.g., block size, levels of fill-in, zero threshold) and (pre-)allocate memory space for $M$ accordingly
  • miscellaneous optimizations (e.g., the Eisenstat trick, see the paper by Chan)
  • code optimizations that target your specific hardware platform (I won't go into specifics here).

At the end of the day, what you're after is a preconditioned iteration that converges faster (in some objective measure) than its un-preconditioned counterpart.

References to get you started:

Yousef Saad, Iterative methods for sparse linear systems - in particular chapters 9.2, 9.2.1 and 9.2.2 deal with the preconditioned CG iteration (there is some pseudocode that you can base your initial implementation on).

Chan and van der Vorst, Approximate and Incomplete Factorizations (they discuss incomplete LU and its variants, but the ideas there apply to ICC as well).

Chih-Jen Lin and Jorge J. Moré, Incomplete Cholesky Factorizations with Limited Memory, SIAM SISC 21(1), 1999

Some high-quality C/C++ ICC implementations (these libraries have a fairly steep learning curve, but the source code is available for inspection):

PETSc's ICC routines.

IFPACK (part of the Trilinos suite). The user manual can be found here.

Eigen's C++ ICC preconditioners for CG.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you for your answer! I just got a few extra question: 1)what do you mean by ‘action vector exactly? I’m afraid that’s pretty unclear to me and 2) what is the block size you mentioned? $\endgroup$
    – Gilles Poncelet
    Commented May 14, 2018 at 13:46
  • 1
    $\begingroup$ (1) by action on a vector I mean the product $M^{-1}v$ for a given vector $v$, that is equivalent to solving the linear system $Mx = v$ for $x$ (again, you do not compute the inverse explicitly). $\endgroup$
    – GoHokies
    Commented May 14, 2018 at 13:52
  • 1
    $\begingroup$ (2) to take advantage of fast cache memory (maximize data reuse), numerical algorithms are often blocked - see here. the optimal block size is dependent on your target architecture. this is an advanced optimization topic that is best ignored at the beginning. $\endgroup$
    – GoHokies
    Commented May 14, 2018 at 14:00
  • $\begingroup$ and I suppose the $Mx = v$ system is super cheap to solve? Which ‘techniques’, if any, are used ? $\endgroup$
    – Gilles Poncelet
    Commented May 14, 2018 at 14:50
  • 1
    $\begingroup$ Well the icc factor $M$ is (upper) triangular (and very sparse) so you can solve that system using backsubstitution $\endgroup$
    – GoHokies
    Commented May 14, 2018 at 14:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.