# Cannot show that the matrix form of my SPD preconditioner is SPD

**URL:** <https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869>\
**Category:** Numerics\
**Created:** [March 10, 2023, 2:53pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869 "2023-03-10T14:53:16Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [March 10, 2023, 2:53pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869/1 "2023-03-10T14:53:16Z")

</div>

I define a symmetric positive definite (SPD) matrix `A`. Then, based on this matrix, I define two preconditioners, which should also be SPD, namely an AMG and a Cholesky factorization. Both preconditioners are matrix-free. Then I define a function `get_matrix_of_precond` which computes the matrix form of each preconditioner. This function is obtained by simply applying each preconditioner to each column of the identity matrix. As a result, I obtain two matrices named `mat_M_amg` and `mat_M_chol`. I expect each of these two matrices to be SPD. However, when I evaluate `isposdef(mat_M_amg)` and `isposdef(mat_M_chol)`, I obtain `false` in both cases. Why is that?

Here is the associated code:

```julia
using LinearAlgebra: cholesky, isposdef
using Preconditioners: AMGPreconditioner, SmoothedAggregation
using SparseArrays: sparse

function get_matrix_of_precond(M, n)
    mat = Array{Float64,2}(undef, n, n)
    for i in 1:n
      vec = zeros(n)
      vec[i] = 1.
      mat[:, i] = M \ vec
    end
    return mat
end

n = 200
X = rand(n, n)
A = X'X

M_amg = AMGPreconditioner{SmoothedAggregation}(sparse(A))
M_chol = cholesky(A)

mat_M_amg = get_matrix_of_precond(M_amg, n)
mat_M_chol = get_matrix_of_precond(M_chol, n)

println("is M_amg posdef? $(isposdef(mat_M_amg))")
println("is M_chol posdef? $(isposdef(mat_M_chol))")

```

---

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [March 10, 2023, 3:09pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869/2 "2023-03-10T15:09:23Z")

</div>

If I check the spectrum of each matrix, they’re both entirely positive, i.e.,

```julia
using LinearAlgebra: eigvals
println(minimum(eigvals(mat_Mat_amg)) > 0)
println(minimum(eigvals(mat_Mat_chol)) > 0)

```

gives me `true` for both cases.

I guess only the symmetry condition is not precisely verified enough when I build the matrix form with `get_matrix_of_precond`. Indeed, if I force the symmetry to be more precisely verified, then `isposdef(mat_M_amg)` and `isposdef(mat_M_chol)` finally return `true`:

```julia
println(isposdef((mat_M_amg' + mat_M_amg) / 2))
println(isposdef((mat_M_chol' + mat_M_chol) / 2))

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 10, 2023, 3:49pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869/3 "2023-03-10T15:49:16Z")

</div>

Another way to ensure symmetry, which brings additional computational advantages, is to store your matrices as a `Symmetric` object

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 10, 2023, 4:00pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869/4 "2023-03-10T16:00:07Z")

</div>

> [@nvenkov1](#):
>
> `function get_matrix_of_precond(M, n)`

PS. Maybe just `get_matrix_of_precond(M) = Symmetric(M \ I(size(M, 1)))`?

---

<div class="post-metadata">

**Author:** ![nvenkov1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nvenkov1/32/43679_2.png) [@nvenkov1](https://discourse.julialang.org/u/nvenkov1)\
**Post date:** [March 10, 2023, 4:12pm UTC](https://discourse.julialang.org/t/cannot-show-that-the-matrix-form-of-my-spd-preconditioner-is-spd/95869/5 "2023-03-10T16:12:50Z")

</div>

`M` can only be applied to vectors. But yeah, I should use `Symmetric`.