I don't think it is a good idea to default to the pivoting version even though it allows the matrix to be positive semidefinite. The pivoting would surprise most user. I don't know your use case but would sqrtm be a possibility instead of chol?
2014-04-08 4:38 GMT+02:00 Jiahao Chen <jia...@mit.edu>: > The behavior of unpivoted Cholesky is not obviously wrong; it fails to > factorize in finite-precision arithmetic when the smallest eigenvalue > is something on the order of epsilon (you can get the precise constant > and citation from, e.g. > > http://mathoverflow.net/questions/108700/quantifying-the-failure-of-the-cholesky-factorization-test-for-indefinite-matric > ). > > I believe the error code returned tells you which row of the > factorization failed. > Thanks, > > Jiahao Chen > Staff Research Scientist > MIT Computer Science and Artificial Intelligence Laboratory > > > On Mon, Apr 7, 2014 at 4:49 PM, Iain Dunning <iaindunn...@gmail.com> > wrote: > > Hey, > > > > Couldn't find the documentation - all I found was > > > https://github.com/JuliaLang/julia/blob/27c3a5b7ed66bee509fc4a81aa54ba09aec1b2ee/base/linalg/exceptions.jl#L19-L21 > > and > > > http://docs.julialang.org/en/latest/stdlib/linalg/?highlight=posdefexception#Base.cholfact > > > > I know rand(4,4) isn't symmetric :) > > I was just trying to figure out if the error code was always 4 or there > was > > some deeper meaning. > > > > My solution is to use cholfact and set pivot to true. Would it be a good > > idea for chol to use pivot by default? Or is the cost unacceptable? > > > > Cheers, > > Iain > > > > > > On Monday, April 7, 2014 4:41:38 PM UTC-4, Douglas Bates wrote: > >> > >> > >> > >> On Monday, April 7, 2014 3:24:32 PM UTC-5, Iain Dunning wrote: > >>> > >>> Hi all, > >>> > >>> I have the following matrix (in a copy-pasteable format) > >>> > >>> X = [1.00000000753845 -.9999999962147147 -1.0000000052345246 > >>> -1.0000000048451771; > >>> -.9999999962147147 .9999999849183356 .9999999938533486 > >>> .9999999932021346; > >>> -1.0000000052345246 .9999999938533486 1.0000000081826723 > >>> 1.0000000034216714; > >>> -1.0000000048451771 .9999999932021346 1.0000000034216714 > >>> 1.0000000058875202] > >>> > >>> and I have > >>> println(eig(X)) > >>> > >>> telling me > >>> [8.01744247656805e-17, 2.6859121865768767e-9, 3.823536365571658e-9, > >>> 4.000000000017529] > >>> (i.e. that it is PSD) > >>> > >>> but > >>> println(chol(X)) > >>> > >>> giving me > >>> ERROR: PosDefException(4) > >>> > >>> Two questions: > >>> 1) How do I make chol try harder? Do I want to use cholfact and > increase > >>> tol to... what? > >>> 2) How can I find out what (4) means? Because when I do the following: > >>> > >> As a last resort you can always read the documentation :-) > >> > >> The 4 is the error number from the LAPACK.potrf! function, indicating > that > >> the Cholesky factorization failed on the fourth column > >> > >>> julia> chol(rand(4,4)) > >>> ERROR: PosDefException(2) > >>> in cholfact! at linalg/factorization.jl:36 > >>> in chol at linalg/factorization.jl:44 > >>> > >> > >> Umm - rand(4,4) is not expected to be symmetric. > >> > >>> > >>> it is (2) - so I'm guessing that the 4 communicates something useful. > >>> > >>> Cheers, > >>> Iain > >>> > >>> > > > -- Med venlig hilsen Andreas Noack Jensen
