2
$\begingroup$

In my case, I am solving $AX=B$ repeatedly, but the solution usually doesn't change much. So it'd probably be faster than me when I start from the previous solution and iterative, rather than solving from scratch. My $A$ is quite dense but it's positive definite.

I see a lot of iterative solver options provided by scipy.. which one should I use in my case?

I am also surprised that they apparently don't support solving matrix $X$ directly, but have to do it by vectors...

https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html

Thanks..

Improve this question
$\endgroup$
10
  • $\begingroup$ What have you tried already? When playing with different options, what do you see? $\endgroup$
    – Wolfgang Bangerth
    Commented Jul 25 at 14:47
  • 1
    $\begingroup$ If the matrix is dense and positive definite and you have many right-hand sides, why not use Cholesky instead? $\endgroup$
    – lightxbulb
    Commented Jul 25 at 16:48
  • $\begingroup$ @lightxbulb if I have a reasonable guess of my solution i'd want to use iterative solver? $\endgroup$
    – Taylor Fang
    Commented Jul 25 at 17:11
  • 1
    $\begingroup$ So is it that you are solving $A_k X_k=B_k$ where the system matrix and right-hand sides change slightly enough so you can reuse the solution $X_k$ as an initial guess for the $k+1$-st step? Then you can just use the conjugate gradient method with the given initial guess. The stopping criterion should then be chosen to not depend on the initial guess and only on $A_k$ and $B_k$. If it is only $B_k$ that changes however, then Cholesky may be more efficient. $\endgroup$
    – lightxbulb
    Commented Jul 25 at 17:17
  • $\begingroup$ What are the dimensions of the matrices? $\endgroup$
    – Royi
    Commented Jul 27 at 7:01

1 Answer 1

Reset to default
2
$\begingroup$

Assuming your model is solving the problem sequentially for slightly different $\boldsymbol{B}_{i}$:

$$ \boldsymbol{A} \boldsymbol{X}_{i} = \boldsymbol{B}_{i} $$

Since the matrix $\boldsymbol{A}$ is a Positive Definite Matrix, I am not sure your approach, using iterative solver, is the most efficient.

I think that computing the Cholesky Decomposition of the matrix and using it for any solution will be faster.

To backup my assertion I created a simple Julia script.
I created 2 matrices of size 5000x5000.

I wanted to compare 2 approaches:

  1. Iterative Approach
    Since your matrix is PD the go to method should be based on Conjugate Gradient. Where the solution of the $i$ -th iterations is the starting point of the $i + 1$ -th solution.
  2. Direct Approach
    Since the matrix is PD, one could cache the Cholesky Decomposoiton and use that for each new right hand matrix.

I neglected the overhead of the first solution of each approach.
Let's focus on the steady state.

In Julia, the (2) approach yields:

runTime = @belapsed ($oChol \ $mB) seconds = 2;
resAnalysis = @sprintf("Direct Solver based on Cholesky run time: %0.5f [Sec]", runTime);
println(resAnalysis);

Direct Solver based on Cholesky run time: 1.04546 [Sec]

Namely, each solution takes something similar to 1 [Sec].

To evaluate the iterative approach, I calculated the time to simply multiply $\boldsymbol{A} \boldsymbol{B}$:

runTime = @belapsed ($mA * $mB) seconds = 2;
resAnalysis = @sprintf("Iterative Solver based on Matrix Multiplication run time: %0.5f [Sec]", runTime);
println(resAnalysis);

Iterative Solver based on Matrix Multiplication run time: 1.06335 [Sec]

So it is on the same ballpark, yet the iterative solver will probably require few iterations and each iteration of CG is using 2 multiplications.

So based on this simple evaluation, it is probably better to stick with the Cholesky based solution.

The code is available on my StackExchange Computational Science GitHub Repository (Look at the ComputationalScience\Q44417 folder).

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ I suspect the winner highly depends on the condition number of the matrix. If CG converges slowly, a good initial value won't make up for it. $\endgroup$
    – Federico Poloni
    Commented Aug 6 at 10:56
  • $\begingroup$ @FedericoPoloni, Since I took random data I'd say the above the is the best case for CG. Don't you think? $\endgroup$
    – Royi
    Commented Aug 6 at 11:07
  • $\begingroup$ I agree that you are in a favorable case for CG, but for a different reason: you generated the matrix with $A = MM^*+I$, so the eigenvalues of $A$ are bounded below by 1, and your matrix is guaranteed to be reasonably well-conditioned (since $||M||$ is moderate). If you change that I to 1e-6 * I I imagine Cholesky will win. (And, on the other hand, if you change it to 1e6 * I I imagine CG will win.) $\endgroup$
    – Federico Poloni
    Commented Aug 6 at 12:01
  • $\begingroup$ @FedericoPoloni, I disagree. Even if CG will converge in 1 iteration (Unreasonable) it would be slower than the Cholesky Decomposition option. Pay attention that my timing compared the operation of Matrix Multiplication. It is still on the same ball park of solving the equation with the decomposition. Hence even a single iteration of CG will be slower path to choose. Even for a single iteration. Hence even for the case the condition number is as close to 1 as possible the CG option will be slower. $\endgroup$
    – Royi
    Commented Aug 6 at 17:54
  • $\begingroup$ Ah, I see what you mean, I missed the fact that mB is a square matrix too. I agree with you, in that setup you can't beat Cholesky. $\endgroup$
    – Federico Poloni
    Commented Aug 6 at 18:12

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.