Re: [julia-users] Re: "Getting around" a PosDefException
Actually much more expensive, but I think that it could easily get confusing to use the permuted Cholesky factor. We should probably have a look at sqrtm and see if it could be made more efficient. It allocates a lot of memory and I am not sure that it is necessary. 2014-04-11 11:28 GMT+02:00 Toivo Henningsson : > Isn't sqrtm more computationally expensive? > > > On Friday, 11 April 2014 09:03:38 UTC+2, Andreas Noack Jensen wrote: >> >> I think that sqrtm would often be the more reasonable advise. My guess is >> that the Cholesky factor is very often used like a matrix square root. >> > -- Med venlig hilsen Andreas Noack Jensen
Re: [julia-users] Re: "Getting around" a PosDefException
Isn't sqrtm more computationally expensive? On Friday, 11 April 2014 09:03:38 UTC+2, Andreas Noack Jensen wrote: > > I think that sqrtm would often be the more reasonable advise. My guess is > that the Cholesky factor is very often used like a matrix square root. >
Re: [julia-users] Re: "Getting around" a PosDefException
I think that sqrtm would often be the more reasonable advise. My guess is that the Cholesky factor is very often used like a matrix square root. 2014-04-08 18:22 GMT+02:00 Stefan Karpinski : > Is it possible to default to unpivoted and if that fails detect that a > pivoted Cholesky might have worked and include a recommendation to try the > pivoted version in the error message? > > > On Tue, Apr 8, 2014 at 10:58 AM, Andreas Noack Jensen < > andreasnoackjen...@gmail.com> wrote: > >> It would be helpful if the LAPACK codes were written out in the Julia >> exception, but it is not most exciting thing to write. The un-pivoted >> Cholesky factor is not triangular, so I think returning that would also >> cause some confusion. >> >> >> 2014-04-08 16:50 GMT+02:00 Iain Dunning : >> >> Jiahao: interesting link! Do you think we should put the meaning of that >>> error code somewhere? Maybe best would be as the actual message of the >>> PosDefException. >>> Andreas: if we un-pivot the result then the user would be unaware, >>> correct? I feel like chol() is the "casual" way of doing it and should make >>> a best effort to work, whereas cholfact is the more poweruser version. >>> David: I was indeed playing around with max-cut, check out >>> https://github.com/JuliaOpt/JuMP.jl/blob/sdp/examples/maxcut_sdp.jl >>> >>> Cheers, >>> Iain >>> >>> >>> On Tuesday, April 8, 2014 5:58:36 AM UTC-4, David de Laat wrote: You can also use a hack to make the matrix positive definite: mineig = minimum(eigvals(M)) M -= mineig * eye(M) (And in case you're working on max-cut you can also use M = (M - mineig * eye(M)) / (1-mineig) so that the linear constraints in the semidefinite program are still satisfied by the new matrix M.) Best, David >>> >> >> >> -- >> Med venlig hilsen >> >> Andreas Noack Jensen >> > > -- Med venlig hilsen Andreas Noack Jensen
Re: [julia-users] Re: "Getting around" a PosDefException
Is it possible to default to unpivoted and if that fails detect that a pivoted Cholesky might have worked and include a recommendation to try the pivoted version in the error message? On Tue, Apr 8, 2014 at 10:58 AM, Andreas Noack Jensen < andreasnoackjen...@gmail.com> wrote: > It would be helpful if the LAPACK codes were written out in the Julia > exception, but it is not most exciting thing to write. The un-pivoted > Cholesky factor is not triangular, so I think returning that would also > cause some confusion. > > > 2014-04-08 16:50 GMT+02:00 Iain Dunning : > > Jiahao: interesting link! Do you think we should put the meaning of that >> error code somewhere? Maybe best would be as the actual message of the >> PosDefException. >> Andreas: if we un-pivot the result then the user would be unaware, >> correct? I feel like chol() is the "casual" way of doing it and should make >> a best effort to work, whereas cholfact is the more poweruser version. >> David: I was indeed playing around with max-cut, check out >> https://github.com/JuliaOpt/JuMP.jl/blob/sdp/examples/maxcut_sdp.jl >> >> Cheers, >> Iain >> >> >> On Tuesday, April 8, 2014 5:58:36 AM UTC-4, David de Laat wrote: >>> >>> You can also use a hack to make the matrix positive definite: >>> mineig = minimum(eigvals(M)) >>> M -= mineig * eye(M) >>> >>> (And in case you're working on max-cut you can also use >>> M = (M - mineig * eye(M)) / (1-mineig) >>> so that the linear constraints in the semidefinite program are still >>> satisfied by the new matrix M.) >>> >>> Best, >>> David >>> >> > > > -- > Med venlig hilsen > > Andreas Noack Jensen >
Re: [julia-users] Re: "Getting around" a PosDefException
It would be helpful if the LAPACK codes were written out in the Julia exception, but it is not most exciting thing to write. The un-pivoted Cholesky factor is not triangular, so I think returning that would also cause some confusion. 2014-04-08 16:50 GMT+02:00 Iain Dunning : > Jiahao: interesting link! Do you think we should put the meaning of that > error code somewhere? Maybe best would be as the actual message of the > PosDefException. > Andreas: if we un-pivot the result then the user would be unaware, > correct? I feel like chol() is the "casual" way of doing it and should make > a best effort to work, whereas cholfact is the more poweruser version. > David: I was indeed playing around with max-cut, check out > https://github.com/JuliaOpt/JuMP.jl/blob/sdp/examples/maxcut_sdp.jl > > Cheers, > Iain > > > On Tuesday, April 8, 2014 5:58:36 AM UTC-4, David de Laat wrote: >> >> You can also use a hack to make the matrix positive definite: >> mineig = minimum(eigvals(M)) >> M -= mineig * eye(M) >> >> (And in case you're working on max-cut you can also use >> M = (M - mineig * eye(M)) / (1-mineig) >> so that the linear constraints in the semidefinite program are still >> satisfied by the new matrix M.) >> >> Best, >> David >> > -- Med venlig hilsen Andreas Noack Jensen
Re: [julia-users] Re: "Getting around" a PosDefException
Jiahao: interesting link! Do you think we should put the meaning of that error code somewhere? Maybe best would be as the actual message of the PosDefException. Andreas: if we un-pivot the result then the user would be unaware, correct? I feel like chol() is the "casual" way of doing it and should make a best effort to work, whereas cholfact is the more poweruser version. David: I was indeed playing around with max-cut, check out https://github.com/JuliaOpt/JuMP.jl/blob/sdp/examples/maxcut_sdp.jl Cheers, Iain On Tuesday, April 8, 2014 5:58:36 AM UTC-4, David de Laat wrote: > > You can also use a hack to make the matrix positive definite: > mineig = minimum(eigvals(M)) > M -= mineig * eye(M) > > (And in case you're working on max-cut you can also use > M = (M - mineig * eye(M)) / (1-mineig) > so that the linear constraints in the semidefinite program are still > satisfied by the new matrix M.) > > Best, > David >
Re: [julia-users] Re: "Getting around" a PosDefException
You can also use a hack to make the matrix positive definite: mineig = minimum(eigvals(M)) M -= mineig * eye(M) (And in case you're working on max-cut you can also use M = (M - mineig * eye(M)) / (1-mineig) so that the linear constraints in the semidefinite program are still satisfied by the new matrix M.) Best, David
Re: [julia-users] Re: "Getting around" a PosDefException
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 : > 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 > 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.753845-.62147147 -1.52345246 > >>> -1.48451771; > >>> -.62147147 .999849183356 .38533486 > >>> .32021346; > >>> -1.52345246 .38533486 1.81826723 > >>> 1.34216714; > >>> -1.48451771 .32021346 1.34216714 > >>> 1.58875202] > >>> > >>> and I have > >>> println(eig(X)) > >>> > >>> telling me > >>> [8.01744247656805e-17, 2.6859121865768767e-9, 3.823536365571658e-9, > >>> 4.0017529] > >>> (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
Re: [julia-users] Re: "Getting around" a PosDefException
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 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.753845-.62147147 -1.52345246 >>> -1.48451771; >>> -.62147147 .999849183356 .38533486 >>> .32021346; >>> -1.52345246 .38533486 1.81826723 >>> 1.34216714; >>> -1.48451771 .32021346 1.34216714 >>> 1.58875202] >>> >>> and I have >>> println(eig(X)) >>> >>> telling me >>> [8.01744247656805e-17, 2.6859121865768767e-9, 3.823536365571658e-9, >>> 4.0017529] >>> (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 >>> >>> >
