# Preconditioner for Block PSD linear system

**URL:** <https://discourse.julialang.org/t/preconditioner-for-block-psd-linear-system/114519>\
**Category:** Numerics\
**Tags:** question, linearalgebra, iterative-solvers\
**Created:** [May 21, 2024, 5:09pm UTC](https://discourse.julialang.org/t/preconditioner-for-block-psd-linear-system/114519 "2024-05-21T17:09:05Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 21, 2024, 5:09pm UTC](https://discourse.julialang.org/t/preconditioner-for-block-psd-linear-system/114519/1 "2024-05-21T17:09:05Z")

</div>

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,

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 21, 2024, 7:09pm UTC](https://discourse.julialang.org/t/preconditioner-for-block-psd-linear-system/114519/2 "2024-05-21T19:09:22Z")

</div>

You can perhaps try the preconditionner P = diag(Sigma\_a^-1, Sigma\_b^1)

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [June 22, 2024, 5:17am UTC](https://discourse.julialang.org/t/preconditioner-for-block-psd-linear-system/114519/3 "2024-06-22T05:17:07Z")

</div>

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:

> **[Reference · LinearOperators.jl](https://jso.dev/LinearOperators.jl/dev/reference/#LinearOperators.BlockDiagonalOperator-Tuple)**
>
> Documentation for 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:

> **[Preconditioners · Krylov.jl](https://jso.dev/Krylov.jl/dev/preconditioners/#Examples)**
>
> Documentation for 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.

```julia
using Krylov

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

```