Struggling with isposdef and floating point errors

General Usage
,
1

I am trying to do a Quasi-Newton method with a basic cholesky adaptation, to guarantee B is positive definite. The code goes as follows:

function choleksyadaption!(A,β = 10e-3, max_iter = 1000)
    @info "Adapting Cholesky Factorization to make Hessian positive definite."
    if (isposdef(A))
        return nothing
    end
    if any(diag(A) .<= 0)
     τ = 0
    else
        τ = minimum(diag(A)) - β
    end

    for i in 1:max_iter
        @info i τ, A
		A = A + τ * I
        if isposdef(A)
            return nothing
        else
            τ = max(2 * τ, β)
        end
    end
    return nothing
end

The problem is, in my eyes, that isposidef, which calls ishermitian requires the non diagonal elements to be exactly the same, which isn’t easy in numerics. For example

ishermitian([1 0.018511830506205723; 0.01851183050620573 1])

returns false, which makes sense, but I am unable to work with it. Is there a workarround, that I am missing?

Thanks!

2

You can use Hermitian to tell the program that a matrix is supposed to be Hermitian. Then it will only look at half the matrix (by default, the upper triangle) and assume that the opposite triangle is conjugate-equal.

julia> using LinearAlgebra

julia> Hermitian([1 2; 3 4], :U) # :U is the default
2×2 Hermitian{Int64, Matrix{Int64}}:
 1  2
 2  4

julia> Hermitian([1 2; 3 4], :L)
2×2 Hermitian{Int64, Matrix{Int64}}:
 1  3
 3  4
3 Likes
3

See also: `isposdef` inconsistent results, numerically unstable? - #10 by andreasnoack