The reason is that two different log methods are used. When doing cholfact!(a) you calculate the Cholesky factorisation in place. The function cholfact! returns a Cholesky type, but the argument remains a Matrix type. Hence when you write logdet(cholfact!(a)) you call the logdet(Cholesky) method whereas cholfact!(a);logged(a) calls the logdet(Matrix) method. The last method doesn't know that the input is in fact the result of a Cholesky factorisation and hence it calculates the logdet of a non-hermitian matrix. The solution is just to bind the result of cholfact! to some variable and then use that variable instead when calling logdet.
Med venlig hilsen Andreas Noack 2014-10-18 19:47 GMT-04:00 Adam Check <[email protected]>: > Hi, > > I am need to compute the cholesky decomposition of a matrix, and then use > the log determinant of the decomposition. I get different answers depending > on if I use cholfact vs. cholfact! > > The code below produces the same number for the logdet (-4.479) > > *srand(1);a=rand(2,2);a=a*a';logdet(cholfact(a))* > > *srand(1);a=rand(2,2);a=a*a';logdet(cholfact!(a))* > > > However, the following code results in different values being reported: > > *srand(1);a=rand(2,2);a=a*a';a=cholfact(a);logdet(a)* > > *srand(1);a=rand(2,2);a=a*a';cholfact!(a);logdet(a)* > with the second code reporting the (incorrect) value of -2.426, even > though the Cholfact object "a" appears identical in each case. > > This error happens consistently when implementing cholfact!() as in the > last line of code, and can result in an error when taking the logdet(): > *ERROR: > DomainError: determinant is negative*, even when no error results from > using the cholfact() function. > > I am currently using julia version 0.3.1, but I believe the problem has > been present in older versions as well. > > - Adam >
