Error by default for pivoted Cholesky of semidefinte matrix

[olof3]

I tried to get the pivoted Cholesky factorization of a positive semidefinite matrix using

C = cholesky(diagm([1, 0]), Val(true))

However this throws RankDeficientException(1), even though it seems like the matrix was factorized alright?!

It took some time for me to figure out that the following works.

C = cholesky(diagm([1, 0]), Val(true). check=false)

To me it seems non-intuitive having to turn off the check for something that should work.

According to the docs one should also be able check issuccess(C), however this throws
ERROR: MethodError: no method matching issuccess(::CholeskyPivoted{Float64,Array{Float64,2}})

Here are the docs:

  cholesky(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

  Compute the pivoted Cholesky factorization of a dense symmetric positive
  semi-definite matrix A and return a CholeskyPivoted factorization. The
  matrix A can either be a Symmetric or Hermitian StridedMatrix or a perfectly
  symmetric or Hermitian StridedMatrix. The triangular Cholesky factor can be
  obtained from the factorization F with: F.L and F.U. The following functions
  are available for CholeskyPivoted objects: size, \, inv, det, and rank. The
  argument tol determines the tolerance for determining the rank. For negative
  values, the tolerance is the machine precision.

  When check = true, an error is thrown if the decomposition fails. When check
  = false, responsibility for checking the decomposition's validity (via
  issuccess) lies with the user.

Another inconsistency is that the docs for LAPACK.pstrf! says that it works for positive-definite matrices while the LAPACK docs says that it works for positive semidefinite matrices. If it's not a typo, a remark in the docs could be in order.

[andreasnoack]

I don’t remember if we did this on purpose but a challenge here is that the info code can be different from zero both when it’s PSD and when it’s indefinite.

I think the issuccess methods should work even though “success” is a bit ambiguous here but that’s the same for LU.

I’d think that our docs for pstrf! are just wrong. Would be great if you could open a PR fixing any or all of these issues.

[olof3]

Hmm, yes, it is a bit confusing that LAPACK gives a non-zero info code for PSD matrices. It seems to have caused some confusion in other places as well.

Does something like

issuccess(F::CholeskyPivoted) = (F.info == 0 || all(F.factors[k,k] >= 0 for k in 1:size(F,1)))

seem reasonable? And using this method to decide if to throw an error if check=true.

[andreasnoack]

Might work but it's been a while since I looked into the detils. The LAPACK documentation says

the matrix A is either rank deficient with computed rank
as returned in RANK, or is not positive semidefinite. See
Section 7 of LAPACK Working Note JuliaLang/julia#161 for further
information.

Would you be able to take a look at section 7?

There is also the consideration of slightly negative diagonal values. We recently introduced a tolerance when computing the matrix square root to avoid that -1e-16 make the square root complex. PSD matrices will probably often end up with a small negative diagonal element.

[olof3]

(here is the working note number 161)

I didn't quite understand if there is an easy way to check that the matrix is PSD. My impression is that there isn't?

"Of course, if there is serious doubt over semi-definiteness then a symmetric indefinite factorization should be used."

In that case I guess the best solution would that the user has to put check=false if s/he known that the matrix is PSD. It should also be made clear from the docs and the error message that this has to be done. And then I guess issuccess(C::CholeskyPivoted) has to check positive definiteness, i.e. info == 0. If it should be intrdouced at all. Not entirely satisfactorily, but the best that can come up with now.