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

Specific Domains Numerics
1

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

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))
3

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

1 Like
4

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

5

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