Two other advantages of just biting the bullet and implementing this for PositiveFactorizations: - there is no iteration step; the first time you try to factor it, it works. - when you think about it, the classic approach of adding a diagonal factor is (in my humble opinion) wrong, wrong, wrong. See the discussion linked from the PositiveFactorization README.
Best, --Tim On Friday, June 10, 2016 5:22:03 PM CDT [email protected] wrote: > Dear Viral and Tim, > > Thanks for your prompt responses to my query. I should have stated more > precisely how I am using positive definiteness testing. The application is > a classic trust-region method for optimization. In a trust region method, > the main operation is as follows. The input is a symmetric (typically > sparse) matrix A and a vector b. The problem is to compute a certain real > parameter lambda. One tests if A + lambda*speye(n) is positive definite; > if so, then one solves a linear system with this coefficient matrix, and if > not, one increases lambda and tries again. > > So the problem with using isposdef is that, in the case that A is actually > positive definite, isposdef discards the Cholesky factor, so then my > application would need to compute it again (redundantly) to solve the > system. In the case of the PostiveFactorizations.jl package, it appears > that the package is aimed at dense rather than sparse matrices. > > So aside from rewriting cholfact, it seems that the only remaining solution > is Viral's suggestion to catch an exception from cholfact. This raises > another question. Most C++ textbooks advise against using exceptions for > ordinary control-flow on the grounds that throwing and catching an > exception is a time-consuming operation. How about in Julia? Is it > reasonable to use try/throw/catch for ordinary control flow in a scientific > code? The Julia manual states that exceptions are much slower than if > statements. But on the other hand, isposdef in cholmod.jl is written in > terms of exceptions! > > Thanks, > Steve > > On Thursday, June 9, 2016 at 10:16:29 PM UTC-4, [email protected] wrote: > > In Matlab to check if a symmetric sparse matrix is positive definite, I > > can say [R,p]=chol(A) and then if p>0 etc. Is this functionality > > available > > in Julia? The cholfact standard routine throws an exception if its > > argument is not positive definite rather than returning any helpful > > information. > > > > I looked at the code for cholfact in cholmod.jl in Base; it appears that I > > can write a modified version of cholfact that exposes this functionality. > > But it would be better if the functionality were available in the library > > so that my code is not sensitive to changes in undocumented low-level > > routines. > > > > Thanks, > > Steve Vavasis
