Logdet of a cholesky! factorized matrix trows a `DomainError`

General Usage
1

The following code trows a DomainError after cholesky!() was called:

mat = [2.4 2.05
       2.05 2.4]

logdet(mat)

logdet(cholesky(mat))

cholesky!(mat)

logdet(mat)

Why is this the case?

2

cholesky! modifies the matrix in place and returns a very thin wrapper excluding the ignored data. The original binding still has that data though:

julia> using LinearAlgebra

julia> mat = [2.4 2.05
              2.05 2.4]
2×2 Array{Float64,2}:
 2.4   2.05
 2.05  2.4

julia> a = cholesky!(mat)
Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 UpperTriangular{Float64,Array{Float64,2}}:
 1.54919  1.32327
  ⋅       0.80558

julia> logdet(a)
0.4430819716794719

julia> mat
2×2 Array{Float64,2}:
 1.54919  1.32327
 2.05     0.80558

Since cholesky! might be called with very large matrices, preserving the existing memory is vital for good performance. Dispatch takes it from there and efficiently works with the wrapper instead of the underlying data.

3

Thanks, so the most performant solution would be (like you did) calling a = cholesky!(mat) and then working with a?

Also, if I am calling isposdef!(mat) to first check whether mat is positive definite, is there a way to further work with the resulting cholesky factorization?

Because in my code, I would like to check whether mat is positive definite and if that is the case, I would like to perform a cholesky factorization.

4

Okay, looking at LinearAlgebra I found that

isposdef!(A::AbstractMatrix) =
    ishermitian(A) && isposdef(cholesky!(Hermitian(A); check = false))

since I know that my matrix is hermitian, I guess the best solution in my case is

mat = [2.4 2.05
       2.05 2.4]

a = cholesky!(Hermitian(mat); check = false)

if isposdef(a)
  ... here I can work with a
else
  ...
1 Like