6
$\begingroup$

I am working on solving Poisson's eq. $\Delta u = -f$ using conjugate gradient method. I am using scipy's linalg.cg function. In this problem, the source function $f$ changes slightly in each iteration, after which I want to solve the equation again. Given that I know that the old solution is close to the current one, how can I use the old solution to create a good preconditioner? That is if $\Delta u_{0} = -f_{0}$ and $\mid f_{0} - f_{1} \mid << 1$, how can I find a good precondition matrix for $\Delta u_{1} = -f_{1}$ for the conjugate gradient method?

$\endgroup$
7
  • 4
    $\begingroup$ Would using $u_0$ as your new starting guess when solving for $u_1$ speed up convergence? $\endgroup$ Apr 19, 2023 at 3:20
  • 1
    $\begingroup$ I'd try using u0 as a deflation space/preconditioner. $\endgroup$ Apr 19, 2023 at 12:55
  • 3
    $\begingroup$ @OmerPaz Do you have problems to store a matrix of size $A$? If not, I would rather precompute + store the LU decomposition. A preconditioner in this case would be non-sense. $\endgroup$
    – ConvexHull
    Apr 19, 2023 at 13:36
  • 3
    $\begingroup$ @OmerPaz I am not an expert, but you may search the topic of "krylov subspace recycling". But I am not quite sure how you could use that with Scipy. $\endgroup$
    – Laurent90
    Apr 20, 2023 at 15:29
  • 1
    $\begingroup$ You can just set the initial guess to the previous solution. Just make sure that the stopping criterion for the implementation of CG that you use is independent of the initial guess so that you may stop earlier. $\endgroup$
    – lightxbulb
    Apr 20, 2023 at 17:37

1 Answer 1

2
$\begingroup$

I will give a short answer considering the LU decomposition instead (see Omer Paz comment)

First of all, we want to solve the following system

$$ \mathbf{A}\vec{x}=\vec{b}. $$

Suppose you are able to store a matrix of size $A$ and your system only changes by (many) different right hand sides $\vec{b}_1,\vec{b}_2, ...,\vec{b}_N$.

A first (naive) idea would be to precompute $\mathbf{A}^{-1}$ and store it, in order to perform (Algorithm 1)

$$ \begin{align} \{\vec{x}_1,\vec{x}_2, ...,\vec{x}_N\}=\mathbf{A}^{-1}\{\vec{b}_1,\vec{b}_2, ...,\vec{b}_N\},\\ \end{align} $$

which is only a matrix vector multiplication.

However, this is not done in practice, since there is an even faster method by precomputing the LU decomposition

$$ \mathbf{L},\mathbf{U},\mathbf{P},\mathbf{Q}\,=\text{lu}\left(\mathbf{A}\right). $$

Here the matrices are defined such that $\mathbf{P}\mathbf{A}\mathbf{Q} = \mathbf{L}\mathbf{U}$, where $\mathbf{L}$ is unit lower triangular, $\mathbf{U}$ is upper triangular, $\mathbf{P}$ is a row permutation matrix/vector and $\mathbf{Q}$ is a column permutation matrix/vector. Now, instead of applying the inverse, you simply perform (Algorithm 2)

$$ \begin{split} \vec{y}~~~~~~&=\text{solve}(\mathbf{L},\,\vec{b}(\mathbf{P})),\\ \vec{x}(\mathbf{Q})&=\text{solve}(\mathbf{U},\,\vec{y}~~~~~), \end{split} $$

which, in most cases, is faster and more accurate than the direct application of $\mathbf{A}^{-1}$ due to memory/cache issues.

The main idea of a preconditioner is to reduce the condition number of a system for iterative solvers with, e.g.

$$ \begin{align} \text{right preconditioner}&&\mathbf{P}&\mathbf{A}&\hspace{-1.7cm}\vec{x}&&=&\hspace{-0.8cm}&\vec{r},\\ \text{left preconditioner}&&&\mathbf{A}\mathbf{P}&\hspace{-1.7cm}\vec{l}&&=&\hspace{-0.8cm}&\vec{b}. \end{align} $$

Now, the best preconditioner you may choose is $\mathbf{P}=\mathbf{A}^{-1}$, since $\mathbf{A}\mathbf{A}^{-1}=\mathbf{A}^{-1}\mathbf{A}=\mathbf{I}$, which would result in an iterative method with only one iteration step (for convergence), since it is equivalent to Algorithm 1. However, as previously stated, there is a even faster Algorithm 2. $\square $

The LU decomposition has to be applied and stored for each system.

$\endgroup$
3
  • 2
    $\begingroup$ I think the question was not on that subject, but rather whether there exist methods to iteratively improve the convergence of krylov methods for consecutive linear system solutions when the matrix is unchanged and it is assumed that the consecutive solutions are close to each other. Even though, of course, LU-decomposion may be an even better solution for sufficiently small systems. $\endgroup$
    – Laurent90
    Apr 20, 2023 at 15:25
  • 3
    $\begingroup$ I agree that this specific answer over-focuses on LU. But I think the sentiment is correct, that if you have an invariant LHS and many RHS, investing in better preconditioners is a great idea and is practically guaranteed to pay dividends (in that it is expected to improve performance regardless of similarity between RHS's). $\endgroup$ Apr 20, 2023 at 15:52
  • $\begingroup$ The answer was motivated by @OmerPaz comment. The answer is not directly related to the original question. $\endgroup$
    – ConvexHull
    Apr 20, 2023 at 16:08

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.