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

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:

``````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))")
``````

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

``````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`:

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

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

1 Like

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

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