Preconditioner for Block PSD linear system

Specific Domains Numerics
, ,
1

Hello,

I am interested in solving the linear system M x = b with an iterative method, where M is positive definite (PD) with the following block structure:

M = \begin{bmatrix} \Sigma_a & 0\\ 0 & \Sigma_b \end{bmatrix} + \begin{bmatrix}A \Sigma_c A^\top & A \Sigma_c B^\top \\ B \Sigma_c A^\top & B \Sigma_c B^\top \end{bmatrix} \\ = \begin{bmatrix} \Sigma_a & 0\\ 0 & \Sigma_b \end{bmatrix} + \begin{bmatrix} A \\ B \end{bmatrix} \Sigma_c \begin{bmatrix}A \\ B \end{bmatrix}^\top

where

  • A and B are sparse linear operators (e.g. discrete differential operators).
  • \Sigma_a and \Sigma_b are PD matrices, such that u \mapsto \Sigma_{a}^{-1} u and v \mapsto \Sigma_{b}^{-1} v are easy to perform.
  • \Sigma_c is only positive semi-definite for which only w \mapsto \Sigma_{c} w can be performed.

I am leaning towards a conjugate gradient iterative solver.
I am wondering if people have experience with this kind of matrix, in particular in the design of preconditioners.

Thank you for your help,

2

You can perhaps try the preconditionner P = diag(Sigma_a^-1, Sigma_b^1)

3

As @rveltz suggested, the block Jacobi preconditioner diag(Sigma_a^{-1}, Sigma_b^{-1}) is a good idea.
I don’t see many alternatives because other blocks are matrix-free operators.

If you already have M as a matrix-free operator, I suggest to use BlockDiagonalOperator from LinearOperators.jl:

You can easily combine your operators to build a preconditioner P = BlockDiagonalOperator(P1, P2).
To build P1 and P2, you can use cholesky from SuiteSparse or LDLFactorization.jl.
I give some examples of them in the documentation of Krylov.jl:

Your operators P1 and P2 can just perform triangular solves with the factors of Sigma_a and Sigma_b.

If.you use Krylov.jl fo the Krylov method, I suggest the following methods: CG, CR, CAR or MINRES.

using Krylov

x, stats = cg(M, b, M=P)
1 Like