Cholesky factorization guarantees on failure?

[rnickb]

If potrf is called on an indefinite matrix A are there any guarantees on what state it leaves A in when it returns?

For example, if the submatrix A'=A[0:k,0:k] is positive definite, but A[0:k+1,0:k+1] is not, can it be assumed that when potrf returns, the A[0:k,0:k] region of memory contains the cholesky factorization of A'?

This looks like it might be true, but I don't see any guarantees listed in the documentation (https://software.intel.com/en-us/node/520881). There are some algorithms that use the partial factorization to compute a bound on the smallest eigenvalue (For example More and Sorensen's algorithm in the paper "Computing a Trust Region Step"), so it would be helpful if such guarantees could be made as part of the contract of potrf

[mecej4]

I'm afraid that, as some of our presidential candidates in the US are saying, the system is rigged :)

A thousand-and-one example runs of indefinite matrices, for which ?POTRF returned a non-zero value for INFO (i.e., Cholesky factorization failed) and yet left the input matrix intact, do not imply a guarantee. A single counter-example, however, should suffice to prove that there can be no guarantee.

Take the SPOTRF example, spotrfx.f, from the MKL distribution, and the matching data file, spotrfx.d. Change the last two diagonal elements to 0.56 and 0.18. You will find that Cholesky factorization fails and the matrix A gets overwritten, even though the leading 3 X 3 submatrix is definite.

The cure, however, is simple, as you probably know: make a back-up copy of A before making the call.

[rnickb]

You missed the point of the question. 

​I know that the potrf function overwrites it's argument matrix. I'm asking what guarantees it makes of the result if the factorization fails because the matrix isn't positive definite. Specifically I'm hoping that it returns with a partially factored (and, I would think, natural form) that allows it to be used for computing the smallest eigenvalue bounds (as in the paper I referenced).

[mecej4]

OK, I see the point now. I do not know if the source code used for the MKL version of ?POTRF is the same as that at Netlib (http://www.netlib.org/lapack/explore-3.1.1-html/dpotrf.f.html), and it is probable that the answer to that is not available to users of MKL.

I looked briefly at the source code at the Netlib link, and the code does not do any pivoting -- as one would expect, pivoting is unnecessary for a definite matrix. On the other hand, it partitions the matrix into smaller blocks and factorizes the blocks. It will take more time to look through the full code (a few hundred lines) and ascertain whether your conjecture (regarding positive definite submatrices being overwritten by their factors even when the full matrix is not positive definite) is correct. 

Do you have the reference for the More-Sorensen paper?

[rnickb]

Here's a link to the More-Sorensen paper http://cds.cern.ch/record/136890/files/CM-P00069193.pdf. The portion I'm referring to is equation (3.9).

Also, described in more detail in this paper by Gay "Computing optimal locally constrained steps" (http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA079719)

[mecej4]

Thanks for the links. I will read the papers in a day or two.

A couple of question, if you don't mind:

   (i) your citing (David M.?) Gay gives a hint that your application may be related to optimization. Is that so?

   (ii) What is the range of matrix sizes of interest? If N <= 128, the Netlib code does an unblocked factorization, which is easier to study.

[rnickb]

Yes, that is correct. The application is solving the subproblem for trust region optimization, which involves repeated cholesky factorizations.

​I'm working with enough variables (>3,000) that the cost of the factorizations are significant, so an optimized factorization routine would be very much beneficial to use, provided it can make the necessary guarantees if the factorization fails.

[Dmitry_B_Intel]

Hi mickb,

Cholesky factorization checks if minor A(1:i,1:i) is positive definite after it has factored minor A(1:i-1,1:i-1).

Documentation says: "If info = i, the leading minor of order i (and therefore the matrix A itself) is not positive-definite, and the factorization could not be completed." This wording does not express the guarantee you requested, namely "If info = i, the upper or lower triangular part of a(1:i-1,1:i-1) is overwritten by the Cholesky factor U or L, as specified by uplo", but the algorithm is constructed so that this property holds.

Thanks
Dima