Determinant and inverse using cholsky parameter

 Classic List Threaded
22 messages
12
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Determinant and inverse using cholsky parameter

 Dear R list members, I have a vector of Cholesky parameterization of a matrix let say A. I would like to compute the determinant and inverse of the original matrix A from the vector of cholesky parameters , would you suggest an R function to do the task. I have tried hard but unable to find anything like that. Please direct me any package or please share R code that can do the task. Thanks and regards, B.Nataraj         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Determinant and inverse using cholsky parameter

 Hi, Isn't the Cholesky decomposition of A=L (L)^T where T stands for "transpose" and L is the Cholesky factor of A. You say you have the  Cholesky decomposition, isn't it L (above)? A<-L%*%t(L) det(A) solve(A) would be your answer. Hope this helps Ozgur
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Determinant and inverse using cholsky parameter

 On Jun 8, 2012, at 16:08 , Özgür Asar wrote: > Hi, > > Isn't the Cholesky decomposition of A=L (L)^T where T stands for "transpose" > and L is the Cholesky factor of A. > > You say you have the  Cholesky decomposition, isn't it L (above)? > > A<-L%*%t(L) > det(A) > solve(A) > > would be your answer. The standard trick is to work directly with L. This is triangular, so the determinant is the product of the diagonal, and the determinant of A is the square of that. Similarly, you can use backsolve to find the inverse of L and multiply it with its transpose for the inverse of A. > > Hope this helps > Ozgur > > > > -- > View this message in context: http://r.789695.n4.nabble.com/Determinant-and-inverse-using-cholsky-parameter-tp4632769p4632808.html> Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Determinant and inverse using cholsky parameter

 In reply to this post by Özgür Asar On Fri, Jun 8, 2012 at 10:08 AM, Özgür Asar <[hidden email]> wrote: > Hi, > > Isn't the Cholesky decomposition of A=L (L)^T where T stands for "transpose" > and L is the Cholesky factor of A. > > You say you have the  Cholesky decomposition, isn't it L (above)? > > A<-L%*%t(L) > det(A) > solve(A) > > would be your answer. A better answer would be to use that L is an (upper or lower, doesn't matter) triangular matrix, so the determinant is just the product of the diagonal elements, so det(A) is just prod(diag(L))^2. kjetil > > Hope this helps > Ozgur > > > > -- > View this message in context: http://r.789695.n4.nabble.com/Determinant-and-inverse-using-cholsky-parameter-tp4632769p4632808.html> Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Determinant and inverse using cholsky parameter

 In reply to this post by nataraj ?chol2inv kjetil On Fri, Jun 8, 2012 at 1:43 AM,  <[hidden email]> wrote: > Dear R list members, > > I have a vector of Cholesky parameterization of a matrix let say A. I would like to compute the determinant and inverse of the original matrix A from the vector of cholesky parameters , would you suggest an R function to do the task. I have tried hard but unable to find anything like that. > > Please direct me any package or please share R code that can do the task. > > Thanks and regards, > B.Nataraj > > >        [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Determinant and inverse using cholsky parameter

 In reply to this post by Özgür Asar Thanks for all who replied me to this topic. I have the "L" (cholesky decomposed matrix)as a vector listed as columwise as mentioned in the article http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.494 and looking for functions to convert the vector into inverse matrix ( but chol2inv mentioned by Kjetil, takes up matrix form) and the corresponding determinant value (square of the product of the diagonal element of the L). I am going to use the function in an maximum likelihood optimization routine for estimating variables in a variance-covariance matrix. Thanks and regards, B.Nataraj -----Original Message----- From: [hidden email] [mailto:[hidden email]] On Behalf Of Özgür Asar Sent: Friday, June 08, 2012 7:38 PM To: [hidden email] Subject: Re: [R] Determinant and inverse using cholsky parameter Hi, Isn't the Cholesky decomposition of A=L (L)^T where T stands for "transpose" and L is the Cholesky factor of A. You say you have the  Cholesky decomposition, isn't it L (above)? A<-L%*%t(L) det(A) solve(A) would be your answer. Hope this helps Ozgur -- View this message in context: http://r.789695.n4.nabble.com/Determinant-and-inverse-using-cholsky-parameter-tp4632769p4632808.htmlSent from the R help mailing list archive at Nabble.com. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Cholesky decomposition error

 In reply to this post by Peter Dalgaard-2 Dear friends, When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error Error in chol.default(M_cov) :   the leading minor of order 10 is not positive definite When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. Could you please help me to find where things are going wrong in my matrix? Thanks and regards, B.Natarj ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Your matrix is not symmetric, positive definite. If you don't know what this means, you shouldn't be using chol() This may be because it isn't to begin with, or due to numerical error, it doesn't behave as one in the decomposition. My relative ignorance of numeric methods for linear algebra prevents me from saying more than that. -- Bert On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: > Dear friends, > > When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error > > Error in chol.default(M_cov) : >  the leading minor of order 10 is not positive definite > > When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. > > Could you please help me to find where things are going wrong in my matrix? > > > Thanks and regards, > B.Natarj > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Hello, Hello, If the input matrix is symmetric, positive definite then the Cholesky decomposition algorithm is stable. That's why it is so used in statistics, where many times the matrices meet those conditions. Therefore, the matrix isn't symmetric, positive definite to begin with. Rui Barradas Em 14-06-2012 13:48, Bert Gunter escreveu: > Your matrix is not symmetric, positive definite. If you don't know > what this means, you shouldn't be using chol() > > This may be because it isn't to begin with, or due to numerical error, > it doesn't behave as one in the decomposition. My relative ignorance > of numeric methods for linear algebra prevents me from saying more > than that. > > -- Bert > > On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >> Dear friends, >> >> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >> >> Error in chol.default(M_cov) : >>   the leading minor of order 10 is not positive definite >> >> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >> >> Could you please help me to find where things are going wrong in my matrix? >> >> >> Thanks and regards, >> B.Natarj >> >> ______________________________________________ >> [hidden email] mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>> and provide commented, minimal, self-contained, reproducible code. > > > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 In reply to this post by Bert Gunter Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned > Error in chol.default(M_cov) : >  the leading minor of order 10 is not positive definite is surfaced during the function call by optim. Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. Regards, B.Nataraj -----Original Message----- From: Bert Gunter [mailto:[hidden email]] Sent: Thursday, June 14, 2012 6:18 PM To: Nataraj B (ORLL-Biotech) Cc: [hidden email] Subject: Re: [R] Cholesky decomposition error Your matrix is not symmetric, positive definite. If you don't know what this means, you shouldn't be using chol() This may be because it isn't to begin with, or due to numerical error, it doesn't behave as one in the decomposition. My relative ignorance of numeric methods for linear algebra prevents me from saying more than that. -- Bert On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: > Dear friends, > > When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error > > Error in chol.default(M_cov) : >  the leading minor of order 10 is not positive definite > > When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. > > Could you please help me to find where things are going wrong in my matrix? > > > Thanks and regards, > B.Natarj > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Follow the posting guide,please: I believe at this point we need reproducible code and your data to provide you help. See ?dput to post your matrix. -- Bert On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: > > Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. > But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned > >> Error in chol.default(M_cov) : >>  the leading minor of order 10 is not positive definite > > is surfaced during the function call by optim. > > Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? > > I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. > > Regards, > B.Nataraj > > > > > > -----Original Message----- > From: Bert Gunter [mailto:[hidden email]] > Sent: Thursday, June 14, 2012 6:18 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > Your matrix is not symmetric, positive definite. If you don't know > what this means, you shouldn't be using chol() > > This may be because it isn't to begin with, or due to numerical error, > it doesn't behave as one in the decomposition. My relative ignorance > of numeric methods for linear algebra prevents me from saying more > than that. > > -- Bert > > On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >> Dear friends, >> >> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >> >> Error in chol.default(M_cov) : >>  the leading minor of order 10 is not positive definite >> >> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >> >> Could you please help me to find where things are going wrong in my matrix? >> >> >> Thanks and regards, >> B.Natarj >> >> ______________________________________________ >> [hidden email] mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>> and provide commented, minimal, self-contained, reproducible code. > > > > -- > > Bert Gunter > Genentech Nonclinical Biostatistics > > Internal Contact Info: > Phone: 467-7374 > Website: > http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm> > -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. This I observed when I change the maximum iteration of the optim function set to 1 and upto iteration no. 3 it runs , it stuck at iteration 4 and above. Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. Regards, B.Nataraj -----Original Message----- From: Bert Gunter [mailto:[hidden email]] Sent: Friday, June 15, 2012 1:51 PM To: Nataraj B (ORLL-Biotech) Cc: [hidden email] Subject: Re: [R] Cholesky decomposition error Follow the posting guide,please: I believe at this point we need reproducible code and your data to provide you help. See ?dput to post your matrix. -- Bert On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: > > Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. > But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned > >> Error in chol.default(M_cov) : >>  the leading minor of order 10 is not positive definite > > is surfaced during the function call by optim. > > Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? > > I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. > > Regards, > B.Nataraj > > > > > > -----Original Message----- > From: Bert Gunter [mailto:[hidden email]] > Sent: Thursday, June 14, 2012 6:18 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > Your matrix is not symmetric, positive definite. If you don't know > what this means, you shouldn't be using chol() > > This may be because it isn't to begin with, or due to numerical error, > it doesn't behave as one in the decomposition. My relative ignorance > of numeric methods for linear algebra prevents me from saying more > than that. > > -- Bert > > On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >> Dear friends, >> >> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >> >> Error in chol.default(M_cov) : >>  the leading minor of order 10 is not positive definite >> >> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >> >> Could you please help me to find where things are going wrong in my matrix? >> >> >> Thanks and regards, >> B.Natarj >> >> ______________________________________________ >> [hidden email] mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>> and provide commented, minimal, self-contained, reproducible code. > > > > -- > > Bert Gunter > Genentech Nonclinical Biostatistics > > Internal Contact Info: > Phone: 467-7374 > Website: > http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm> > -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 see inline. On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: > Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. ¿What are you trying to optimize? If the objective function depends on a matrix argument which has to be a positive definite function, you must parametrize the matrix such that the matrix inside the optimizer always is positive definite. So if your positive definite matrix is A, then, for example, represent it as its cholesky decomposition A= L L^T where L is lower triangular with positive diagonal. Here the stricly upper diagonal part varies freely, but the diagonal not, so represent the diagonal as exp( l_i) where now the l_i varies freely. This is called the log-Cholesky parametrization. For other ideas along this lines, see the paper by Douglas Bates:  "Unconstrained Parameterizations for Variance-Covariance Matrices  " which you can find by googling. Kjetil  This I observed when I change the maximum iteration of the optim function set to 1 and upto iteration no. 3 it runs , it stuck at iteration 4 and above. > > Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. > > Regards, > B.Nataraj > > > -----Original Message----- > From: Bert Gunter [mailto:[hidden email]] > Sent: Friday, June 15, 2012 1:51 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > Follow the posting guide,please: I believe at this point we need > reproducible code and your data to provide you help. See ?dput to post > your matrix. > > -- Bert > > On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >> >> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >> >>> Error in chol.default(M_cov) : >>>  the leading minor of order 10 is not positive definite >> >> is surfaced during the function call by optim. >> >> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >> >> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >> >> Regards, >> B.Nataraj >> >> >> >> >> >> -----Original Message----- >> From: Bert Gunter [mailto:[hidden email]] >> Sent: Thursday, June 14, 2012 6:18 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> Your matrix is not symmetric, positive definite. If you don't know >> what this means, you shouldn't be using chol() >> >> This may be because it isn't to begin with, or due to numerical error, >> it doesn't behave as one in the decomposition. My relative ignorance >> of numeric methods for linear algebra prevents me from saying more >> than that. >> >> -- Bert >> >> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>> Dear friends, >>> >>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>> >>> Error in chol.default(M_cov) : >>>  the leading minor of order 10 is not positive definite >>> >>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>> >>> Could you please help me to find where things are going wrong in my matrix? >>> >>> >>> Thanks and regards, >>> B.Natarj >>> >>> ______________________________________________ >>> [hidden email] mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>> and provide commented, minimal, self-contained, reproducible code. >> >> >> >> -- >> >> Bert Gunter >> Genentech Nonclinical Biostatistics >> >> Internal Contact Info: >> Phone: 467-7374 >> Website: >> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>> >> > > > > -- > > Bert Gunter > Genentech Nonclinical Biostatistics > > Internal Contact Info: > Phone: 467-7374 > Website: > http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm> > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Dear Mr.Kjetil, Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. Regards, B.Nataraj -----Original Message----- From: Kjetil Halvorsen [mailto:[hidden email]] Sent: Friday, June 15, 2012 11:09 PM To: Nataraj B (ORLL-Biotech) Cc: [hidden email]; [hidden email] Subject: Re: [R] Cholesky decomposition error see inline. On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: > Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. ¿What are you trying to optimize? If the objective function depends on a matrix argument which has to be a positive definite function, you must parametrize the matrix such that the matrix inside the optimizer always is positive definite. So if your positive definite matrix is A, then, for example, represent it as its cholesky decomposition A= L L^T where L is lower triangular with positive diagonal. Here the stricly upper diagonal part varies freely, but the diagonal not, so represent the diagonal as exp( l_i) where now the l_i varies freely. This is called the log-Cholesky parametrization. For other ideas along this lines, see the paper by Douglas Bates:  "Unconstrained Parameterizations for Variance-Covariance Matrices  " which you can find by googling. Kjetil  This I observed when I change the maximum iteration of the optim function set to 1 and upto iteration no. 3 it runs , it stuck at iteration 4 and above. > > Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. > > Regards, > B.Nataraj > > > -----Original Message----- > From: Bert Gunter [mailto:[hidden email]] > Sent: Friday, June 15, 2012 1:51 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > Follow the posting guide,please: I believe at this point we need > reproducible code and your data to provide you help. See ?dput to post > your matrix. > > -- Bert > > On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >> >> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >> >>> Error in chol.default(M_cov) : >>>  the leading minor of order 10 is not positive definite >> >> is surfaced during the function call by optim. >> >> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >> >> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >> >> Regards, >> B.Nataraj >> >> >> >> >> >> -----Original Message----- >> From: Bert Gunter [mailto:[hidden email]] >> Sent: Thursday, June 14, 2012 6:18 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> Your matrix is not symmetric, positive definite. If you don't know >> what this means, you shouldn't be using chol() >> >> This may be because it isn't to begin with, or due to numerical error, >> it doesn't behave as one in the decomposition. My relative ignorance >> of numeric methods for linear algebra prevents me from saying more >> than that. >> >> -- Bert >> >> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>> Dear friends, >>> >>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>> >>> Error in chol.default(M_cov) : >>>  the leading minor of order 10 is not positive definite >>> >>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>> >>> Could you please help me to find where things are going wrong in my matrix? >>> >>> >>> Thanks and regards, >>> B.Natarj >>> >>> ______________________________________________ >>> [hidden email] mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>> and provide commented, minimal, self-contained, reproducible code. >> >> >> >> -- >> >> Bert Gunter >> Genentech Nonclinical Biostatistics >> >> Internal Contact Info: >> Phone: 467-7374 >> Website: >> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>> >> > > > > -- > > Bert Gunter > Genentech Nonclinical Biostatistics > > Internal Contact Info: > Phone: 467-7374 > Website: > http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm> > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 see below. On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: > Dear Mr.Kjetil, > > Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. > > I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. > > Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. I dont think there is a prewritten function for log cholesy, but it is very easy to write! > makelogchol function(A) {#A should be a square posdef matrix m <- dim(A)[1] uind <- upper.tri(A) R <- chol(A) diag <- log(diag(R)) upper <- R[uind] return(list(diag=diag, upper=upper)) } and its "inverse" --- to get the cholesky factor back: > makemat function(diag, upper) { m <- length(diag) mu <- m*(m-1)/2 if (length(upper)!=mu) stop("incompatible lengths") A <- matrix(0.0, m, m) ind <- upper.tri(A) A[ind] <- upper diag(A) <- exp(diag) A } and an example of its use: Lets say A is the posdef matrix where you will start your optimization, so you need the log cholesy parameters to feed to optim: > A      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > start <- makelogchol(A) > start \$diag [1] 0.5493061 0.4904146 0.4581454 \$upper [1] 0.5773503 0.5773503 0.4082483 Then to get the cholesky factors back you call makemat: > R <- makemat(start\$diag,  start\$upper) > R          [,1]      [,2]      [,3] [1,] 1.732051 0.5773503 0.5773503 [2,] 0.000000 1.6329932 0.4082483 [3,] 0.000000 0.0000000 1.5811388 > t(R) %*% R      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > Kjetil > > Regards, > B.Nataraj > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Friday, June 15, 2012 11:09 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email]; [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see inline. > > On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. > > ¿What are you trying to optimize? If the objective function depends on > a matrix argument > which has to be a positive definite function, you must parametrize the > matrix such that the matrix > inside the optimizer always is positive definite. So if your positive > definite matrix is A, then, for example, represent it as > its cholesky decomposition A= L L^T where L is lower triangular with > positive diagonal. Here the stricly > upper diagonal part varies freely, but the diagonal not, so represent > the diagonal as exp( l_i) > where now the l_i varies freely. This is called the log-Cholesky > parametrization. For other ideas along this lines, see the paper by > Douglas Bates:  "Unconstrained Parameterizations for > Variance-Covariance Matrices  " > which you can find by googling. > > Kjetil > > >  This I observed when I change the maximum iteration of the optim > function set to 1 and upto iteration no. 3 it runs , it stuck at > iteration 4 and above. >> >> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >> >> Regards, >> B.Nataraj >> >> >> -----Original Message----- >> From: Bert Gunter [mailto:[hidden email]] >> Sent: Friday, June 15, 2012 1:51 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> Follow the posting guide,please: I believe at this point we need >> reproducible code and your data to provide you help. See ?dput to post >> your matrix. >> >> -- Bert >> >> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>> >>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>> >>>> Error in chol.default(M_cov) : >>>>  the leading minor of order 10 is not positive definite >>> >>> is surfaced during the function call by optim. >>> >>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>> >>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>> >>> Regards, >>> B.Nataraj >>> >>> >>> >>> >>> >>> -----Original Message----- >>> From: Bert Gunter [mailto:[hidden email]] >>> Sent: Thursday, June 14, 2012 6:18 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> Your matrix is not symmetric, positive definite. If you don't know >>> what this means, you shouldn't be using chol() >>> >>> This may be because it isn't to begin with, or due to numerical error, >>> it doesn't behave as one in the decomposition. My relative ignorance >>> of numeric methods for linear algebra prevents me from saying more >>> than that. >>> >>> -- Bert >>> >>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>> Dear friends, >>>> >>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>> >>>> Error in chol.default(M_cov) : >>>>  the leading minor of order 10 is not positive definite >>>> >>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>> >>>> Could you please help me to find where things are going wrong in my matrix? >>>> >>>> >>>> Thanks and regards, >>>> B.Natarj >>>> >>>> ______________________________________________ >>>> [hidden email] mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> >>> >>> -- >>> >>> Bert Gunter >>> Genentech Nonclinical Biostatistics >>> >>> Internal Contact Info: >>> Phone: 467-7374 >>> Website: >>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>> >>> >> >> >> >> -- >> >> Bert Gunter >> Genentech Nonclinical Biostatistics >> >> Internal Contact Info: >> Phone: 467-7374 >> Website: >> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>> >> ______________________________________________ >> [hidden email] mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>> and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Thanks Kjetil for your detailed code for log-Cholesky but my situation is like optimization of the variables inside the matrix. Let me change your own example matrix to explain the problem.  A martrix is      [,1]         [,2]            [,3]                 [1,]    3*Vs+Vn    1*Vs         1*Vs                 [2,]    1*Vs     3*Vs+Vn    1*Vs                 [3,]    1*Vs     1*Vs       3*Vs+Vn So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood function, which need the determinant and inverse of the matrix A to compute the value. Isn't it logically correct that I have to seed some initial values for the Vs and Vn and use it for predicting the value of the matrix and check for its positive definite condition, if the condition fulfill, then I have to compute the determinant value and inverse of the matrix in order to use them in the maximum likelihood function and the optimization iteration to be carried out further until parameters converges. I sincerely believe that you have points to implements log-Cholesky in my situation, but I am simply not able to see where is the Choleksy decomposition to be implemented in this my workflow. Regards, B.Nataraj -----Original Message----- From: Kjetil Halvorsen [mailto:[hidden email]] Sent: Sunday, June 17, 2012 4:10 AM To: Nataraj B (ORLL-Biotech) Cc: [hidden email] Subject: Re: [R] Cholesky decomposition error see below. On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: > Dear Mr.Kjetil, > > Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. > > I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. > > Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. I dont think there is a prewritten function for log cholesy, but it is very easy to write! > makelogchol function(A) {#A should be a square posdef matrix m <- dim(A)[1] uind <- upper.tri(A) R <- chol(A) diag <- log(diag(R)) upper <- R[uind] return(list(diag=diag, upper=upper)) } and its "inverse" --- to get the cholesky factor back: > makemat function(diag, upper) { m <- length(diag) mu <- m*(m-1)/2 if (length(upper)!=mu) stop("incompatible lengths") A <- matrix(0.0, m, m) ind <- upper.tri(A) A[ind] <- upper diag(A) <- exp(diag) A } and an example of its use: Lets say A is the posdef matrix where you will start your optimization, so you need the log cholesy parameters to feed to optim: > A      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > start <- makelogchol(A) > start \$diag [1] 0.5493061 0.4904146 0.4581454 \$upper [1] 0.5773503 0.5773503 0.4082483 Then to get the cholesky factors back you call makemat: > R <- makemat(start\$diag,  start\$upper) > R          [,1]      [,2]      [,3] [1,] 1.732051 0.5773503 0.5773503 [2,] 0.000000 1.6329932 0.4082483 [3,] 0.000000 0.0000000 1.5811388 > t(R) %*% R      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > Kjetil > > Regards, > B.Nataraj > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Friday, June 15, 2012 11:09 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email]; [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see inline. > > On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. > > ¿What are you trying to optimize? If the objective function depends on > a matrix argument > which has to be a positive definite function, you must parametrize the > matrix such that the matrix > inside the optimizer always is positive definite. So if your positive > definite matrix is A, then, for example, represent it as > its cholesky decomposition A= L L^T where L is lower triangular with > positive diagonal. Here the stricly > upper diagonal part varies freely, but the diagonal not, so represent > the diagonal as exp( l_i) > where now the l_i varies freely. This is called the log-Cholesky > parametrization. For other ideas along this lines, see the paper by > Douglas Bates:  "Unconstrained Parameterizations for > Variance-Covariance Matrices  " > which you can find by googling. > > Kjetil > > >  This I observed when I change the maximum iteration of the optim > function set to 1 and upto iteration no. 3 it runs , it stuck at > iteration 4 and above. >> >> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >> >> Regards, >> B.Nataraj >> >> >> -----Original Message----- >> From: Bert Gunter [mailto:[hidden email]] >> Sent: Friday, June 15, 2012 1:51 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> Follow the posting guide,please: I believe at this point we need >> reproducible code and your data to provide you help. See ?dput to post >> your matrix. >> >> -- Bert >> >> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>> >>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>> >>>> Error in chol.default(M_cov) : >>>>  the leading minor of order 10 is not positive definite >>> >>> is surfaced during the function call by optim. >>> >>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>> >>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>> >>> Regards, >>> B.Nataraj >>> >>> >>> >>> >>> >>> -----Original Message----- >>> From: Bert Gunter [mailto:[hidden email]] >>> Sent: Thursday, June 14, 2012 6:18 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> Your matrix is not symmetric, positive definite. If you don't know >>> what this means, you shouldn't be using chol() >>> >>> This may be because it isn't to begin with, or due to numerical error, >>> it doesn't behave as one in the decomposition. My relative ignorance >>> of numeric methods for linear algebra prevents me from saying more >>> than that. >>> >>> -- Bert >>> >>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>> Dear friends, >>>> >>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>> >>>> Error in chol.default(M_cov) : >>>>  the leading minor of order 10 is not positive definite >>>> >>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>> >>>> Could you please help me to find where things are going wrong in my matrix? >>>> >>>> >>>> Thanks and regards, >>>> B.Natarj >>>> >>>> ______________________________________________ >>>> [hidden email] mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> >>> >>> -- >>> >>> Bert Gunter >>> Genentech Nonclinical Biostatistics >>> >>> Internal Contact Info: >>> Phone: 467-7374 >>> Website: >>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>> >>> >> >> >> >> -- >> >> Bert Gunter >> Genentech Nonclinical Biostatistics >> >> Internal Contact Info: >> Phone: 467-7374 >> Website: >> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>> >> ______________________________________________ >> [hidden email] mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>> and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 see inline! On Mon, Jun 18, 2012 at 12:38 AM,  <[hidden email]> wrote: > Thanks Kjetil for your detailed code for log-Cholesky but my situation is like optimization of the variables inside the matrix. But you must of course call my function makemat inside the call to optim() or whatever function you use for optimization! Here is an extended example, assuming the functions I have defined earlier: library(MASS) # for mvrnorm library(mvtnorm) # for dmvnorm Let us simulate a multinorma sample with expectation zero (to simplify the example) and covariance matrix Sigma: > Sigma      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > > testdat <- mvrnorm(n=100, rep(0, 3), Sigma) We can calculate the usual maxlik estimator of Sigma: > (t(testdat) %*% testdat) /100           [,1]      [,2]     [,3] [1,] 3.2424618 0.9670133 1.673897 [2,] 0.9670133 2.9447654 1.289127 [3,] 1.6738970 1.2891268 3.084106 > objective <- function(Sigma) {         sum(apply(testdat, 1, FUN=function(x) dmvnorm(x, rep(0, 3), Sigma, log=TRUE)))     } > objective(Sigma) > start \$diag [1] 0.5493061 0.4904146 0.4581454 \$upper [1] 0.5773503 0.5773503 0.4082483 > startSigma <- c(start\$diag, start\$upper) > startSigma [1] 0.5493061 0.4904146 0.4581454 0.5773503 0.5773503 0.4082483 > objective(Sigma) [1] -571.5945 > objective(Sigma/3) [1] -699.0552 > objective(3*Sigma) [1] -638.9688 > Finally:  optim(startSigma, fn=function(parvec) {         R <- makemat(parvec[1:3], parvec[4:6])         S <- t(R) %*% R         print(S)   # remove this when works!        obj <- -objective(S)        obj     }, control=list(REPORT=1)) # with method="BFGS" dopes not work well,  because takes too long # steps! default Nelder-Mead is slow,  but works. Giving the following output: <  much removed ...> \$par [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 \$value [1] 567.648 \$counts function gradient      501       NA \$convergence [1] 1 \$message NULL > > res <- .Last.value > res\$par [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 > makemat(res\$par[1:3], res\$par[4:6])          [,1]     [,2]      [,3] [1,] 1.788146 0.546956 0.9248746 [2,] 0.000000 1.643628 0.4896758 [3,] 0.000000 0.000000 1.4170710 > R <-makemat(res\$par[1:3], res\$par[4:6]) > t(R) %*% R           [,1]      [,2]     [,3] [1,] 3.1974676 0.9780374 1.653811 [2,] 0.9780374 3.0006752 1.310711 [3,] 1.6538112 1.3107108 3.103266 > Sigma      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 ¿Easy? > > Let me change your own example matrix to explain the problem. > >  A martrix is      [,1]         [,2]            [,3] >                [1,]    3*Vs+Vn    1*Vs         1*Vs >                [2,]    1*Vs     3*Vs+Vn    1*Vs >                [3,]    1*Vs     1*Vs       3*Vs+Vn So make your own makemat function specific to this problem! which you then use inside the call to optim. kjetil > > So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood function, which need the determinant and inverse of the matrix A to compute the value. > > Isn't it logically correct that I have to seed some initial values for the Vs and Vn and use it for predicting the value of the matrix and check for its positive definite condition, If you construct the matrix by above method, it will automatically be posdef, so no need to check! (except as check on programming ... ) if the condition fulfill, then I have to compute the determinant value and inverse of the matrix in order to use them in the maximum likelihood function and the optimization iteration to be carried out further until parameters converges. > > I sincerely believe that you have points to implements log-Cholesky in my situation, but I am simply not able to see where is the Choleksy decomposition to be implemented in this my workflow. > > Regards, > B.Nataraj > > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Sunday, June 17, 2012 4:10 AM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see below. > > On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: >> Dear Mr.Kjetil, >> >> Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. >> >> I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. >> >> Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. > > I dont think there is a prewritten function for log cholesy, but it is > very easy to write! > >> makelogchol > function(A) {#A should be a square posdef matrix > m <- dim(A)[1] > uind <- upper.tri(A) > R <- chol(A) > diag <- log(diag(R)) > upper <- R[uind] > return(list(diag=diag, upper=upper)) > } > > and its "inverse" --- to get the cholesky factor back: > >> makemat > function(diag, upper) { > m <- length(diag) > mu <- m*(m-1)/2 > if (length(upper)!=mu) stop("incompatible lengths") > A <- matrix(0.0, m, m) > ind <- upper.tri(A) > A[ind] <- upper > diag(A) <- exp(diag) > A > } > > and an example of its use: > > Lets say A is the posdef matrix where you will start your > optimization, so you need > the log cholesy parameters to feed to optim: > >> A >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> start <- makelogchol(A) >> start > \$diag > [1] 0.5493061 0.4904146 0.4581454 > > \$upper > [1] 0.5773503 0.5773503 0.4082483 > > Then to get the cholesky factors back you call makemat: > >> R <- makemat(start\$diag,  start\$upper) >> R >         [,1]      [,2]      [,3] > [1,] 1.732051 0.5773503 0.5773503 > [2,] 0.000000 1.6329932 0.4082483 > [3,] 0.000000 0.0000000 1.5811388 >> t(R) %*% R >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> > > Kjetil > > > > > > > >> >> Regards, >> B.Nataraj >> >> >> >> -----Original Message----- >> From: Kjetil Halvorsen [mailto:[hidden email]] >> Sent: Friday, June 15, 2012 11:09 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email]; [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> see inline. >> >> On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >>> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. >> >> ¿What are you trying to optimize? If the objective function depends on >> a matrix argument >> which has to be a positive definite function, you must parametrize the >> matrix such that the matrix >> inside the optimizer always is positive definite. So if your positive >> definite matrix is A, then, for example, represent it as >> its cholesky decomposition A= L L^T where L is lower triangular with >> positive diagonal. Here the stricly >> upper diagonal part varies freely, but the diagonal not, so represent >> the diagonal as exp( l_i) >> where now the l_i varies freely. This is called the log-Cholesky >> parametrization. For other ideas along this lines, see the paper by >> Douglas Bates:  "Unconstrained Parameterizations for >> Variance-Covariance Matrices  " >> which you can find by googling. >> >> Kjetil >> >> >>  This I observed when I change the maximum iteration of the optim >> function set to 1 and upto iteration no. 3 it runs , it stuck at >> iteration 4 and above. >>> >>> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >>> >>> Regards, >>> B.Nataraj >>> >>> >>> -----Original Message----- >>> From: Bert Gunter [mailto:[hidden email]] >>> Sent: Friday, June 15, 2012 1:51 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> Follow the posting guide,please: I believe at this point we need >>> reproducible code and your data to provide you help. See ?dput to post >>> your matrix. >>> >>> -- Bert >>> >>> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>>> >>>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>>> >>>>> Error in chol.default(M_cov) : >>>>>  the leading minor of order 10 is not positive definite >>>> >>>> is surfaced during the function call by optim. >>>> >>>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>>> >>>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>>> >>>> Regards, >>>> B.Nataraj >>>> >>>> >>>> >>>> >>>> >>>> -----Original Message----- >>>> From: Bert Gunter [mailto:[hidden email]] >>>> Sent: Thursday, June 14, 2012 6:18 PM >>>> To: Nataraj B (ORLL-Biotech) >>>> Cc: [hidden email] >>>> Subject: Re: [R] Cholesky decomposition error >>>> >>>> Your matrix is not symmetric, positive definite. If you don't know >>>> what this means, you shouldn't be using chol() >>>> >>>> This may be because it isn't to begin with, or due to numerical error, >>>> it doesn't behave as one in the decomposition. My relative ignorance >>>> of numeric methods for linear algebra prevents me from saying more >>>> than that. >>>> >>>> -- Bert >>>> >>>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>>> Dear friends, >>>>> >>>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>>> >>>>> Error in chol.default(M_cov) : >>>>>  the leading minor of order 10 is not positive definite >>>>> >>>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>>> >>>>> Could you please help me to find where things are going wrong in my matrix? >>>>> >>>>> >>>>> Thanks and regards, >>>>> B.Natarj >>>>> >>>>> ______________________________________________ >>>>> [hidden email] mailing list >>>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>>> and provide commented, minimal, self-contained, reproducible code. >>>> >>>> >>>> >>>> -- >>>> >>>> Bert Gunter >>>> Genentech Nonclinical Biostatistics >>>> >>>> Internal Contact Info: >>>> Phone: 467-7374 >>>> Website: >>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>> >>>> >>> >>> >>> >>> -- >>> >>> Bert Gunter >>> Genentech Nonclinical Biostatistics >>> >>> Internal Contact Info: >>> Phone: 467-7374 >>> Website: >>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>> >>> ______________________________________________ >>> [hidden email] mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>> and provide commented, minimal, self-contained, reproducible code. >> >> > > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Dear Kjetil, thanks for your time to detail the code. Sorry, I am still not able to grasp why you are creating random normal distribution for the covariance matrix and also why I should optimize the whole matrix value while I only need to optimize the two variables ?. This is my 15 x 15 covariance matrix > M_cov [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]  [1,]    4    3    3    2    2    1    2    2    3     2     2     2     3     3     1  [2,]    3    5    3    2    2    1    2    2    2     2     2     2     2     2     1  [3,]    3    3    3    2    2    1    2    2    2     2     2     2     2     2     1  [4,]    2    2    2    3    2    1    2    3    2     1     1     1     1     1     1  [5,]    2    2    2    2    2    1    2    2    2     1     1     1     1     1     1  [6,]    1    1    1    1    1    1    1    1    1     1     0     1     1     0     0  [7,]    2    2    2    2    2    1    3    2    2     1     1     1     1     1     1  [8,]    2    2    2    3    2    1    2    4    2     1     1     1     1     1     1  [9,]    3    2    2    2    2    1    2    2    3     1     1     1     2     2     1 [10,]    2    2    2    1    1    1    1    1    1     2     1     2     2     1     0 [11,]    2    2    2    1    1    0    1    1    1     1     2     1     1     2     1 [12,]    2    2    2    1    1    1    1    1    1     2     1     3     2     1     0 [13,]    3    2    2    1    1    1    1    1    2     2     1     2     3     2     0 [14,]    3    2    2    1    1    0    1    1    2     1     2     1     2     3     1 [15,]    1    1    1    1    1    0    1    1    1     0     1     0     0     1     1 And the function to check the posdef condition is >isPosDef <- function(M) { if ( all(M == t(M) ) ) {  # first test symmetric-ity   if (  all(eigen(M)\$values > 0) ) {TRUE}     else {FALSE} } #  else {FALSE}  # not symmetric  } And my call for optimization is >optim(c(rep(1,2)),obj,method="BFGS",hessian=TRUE,F1=F1,mu1=mu1,M_cov=M_cov) And I pass the vector F1 and mu1 used in my objective function along with M_cov (the covariance matrix) My objective function is >obj<-function(theta,F1,mu1,M_cov) { Vs<-theta[1] Vn<-theta[2] k<-length(F1) for(i in 1:15) { for(j in 1:15) { # to make the covariance matrix with the variables to optimize if(i != j) { M_cov[i,j]<-M_cov[i,j]*Vs } else M_cov[i,j]<-M_cov[i,j]*Vs+Vn } } # here I check the posdef condition of the modified covariance matrix in #earlier step checkPD<-isPosDef(M_cov) if(checkPD) { e<-t(F1-mu1)%*%solve(M_cov)%*%(F1-mu1) logl<- -k/2*log(2*pi)-0.5*log(det(M_cov))-0.5*e } # negative value since the function is for maximization return(-logl) } Is there anything missing in my approach of function optimization for variable estimation of a multivariate normal distribution? I looking forward your valuable comments. Regards, B.Nataraj -----Original Message----- From: Kjetil Halvorsen [mailto:[hidden email]] Sent: Monday, June 18, 2012 10:07 PM To: Nataraj B (ORLL-Biotech) Cc: [hidden email] Subject: Re: [R] Cholesky decomposition error see inline! On Mon, Jun 18, 2012 at 12:38 AM,  <[hidden email]> wrote: > Thanks Kjetil for your detailed code for log-Cholesky but my situation is like optimization of the variables inside the matrix. But you must of course call my function makemat inside the call to optim() or whatever function you use for optimization! Here is an extended example, assuming the functions I have defined earlier: library(MASS) # for mvrnorm library(mvtnorm) # for dmvnorm Let us simulate a multinorma sample with expectation zero (to simplify the example) and covariance matrix Sigma: > Sigma      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 > > testdat <- mvrnorm(n=100, rep(0, 3), Sigma) We can calculate the usual maxlik estimator of Sigma: > (t(testdat) %*% testdat) /100           [,1]      [,2]     [,3] [1,] 3.2424618 0.9670133 1.673897 [2,] 0.9670133 2.9447654 1.289127 [3,] 1.6738970 1.2891268 3.084106 > objective <- function(Sigma) {         sum(apply(testdat, 1, FUN=function(x) dmvnorm(x, rep(0, 3), Sigma, log=TRUE)))     } > objective(Sigma) > start \$diag [1] 0.5493061 0.4904146 0.4581454 \$upper [1] 0.5773503 0.5773503 0.4082483 > startSigma <- c(start\$diag, start\$upper) > startSigma [1] 0.5493061 0.4904146 0.4581454 0.5773503 0.5773503 0.4082483 > objective(Sigma) [1] -571.5945 > objective(Sigma/3) [1] -699.0552 > objective(3*Sigma) [1] -638.9688 > Finally:  optim(startSigma, fn=function(parvec) {         R <- makemat(parvec[1:3], parvec[4:6])         S <- t(R) %*% R         print(S)   # remove this when works!        obj <- -objective(S)        obj     }, control=list(REPORT=1)) # with method="BFGS" dopes not work well,  because takes too long # steps! default Nelder-Mead is slow,  but works. Giving the following output: <  much removed ...> \$par [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 \$value [1] 567.648 \$counts function gradient      501       NA \$convergence [1] 1 \$message NULL > > res <- .Last.value > res\$par [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 > makemat(res\$par[1:3], res\$par[4:6])          [,1]     [,2]      [,3] [1,] 1.788146 0.546956 0.9248746 [2,] 0.000000 1.643628 0.4896758 [3,] 0.000000 0.000000 1.4170710 > R <-makemat(res\$par[1:3], res\$par[4:6]) > t(R) %*% R           [,1]      [,2]     [,3] [1,] 3.1974676 0.9780374 1.653811 [2,] 0.9780374 3.0006752 1.310711 [3,] 1.6538112 1.3107108 3.103266 > Sigma      [,1] [,2] [,3] [1,]    3    1    1 [2,]    1    3    1 [3,]    1    1    3 ¿Easy? > > Let me change your own example matrix to explain the problem. > >  A martrix is      [,1]         [,2]            [,3] >                [1,]    3*Vs+Vn    1*Vs         1*Vs >                [2,]    1*Vs     3*Vs+Vn    1*Vs >                [3,]    1*Vs     1*Vs       3*Vs+Vn So make your own makemat function specific to this problem! which you then use inside the call to optim. kjetil > > So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood function, which need the determinant and inverse of the matrix A to compute the value. > > Isn't it logically correct that I have to seed some initial values for the Vs and Vn and use it for predicting the value of the matrix and check for its positive definite condition, If you construct the matrix by above method, it will automatically be posdef, so no need to check! (except as check on programming ... ) if the condition fulfill, then I have to compute the determinant value and inverse of the matrix in order to use them in the maximum likelihood function and the optimization iteration to be carried out further until parameters converges. > > I sincerely believe that you have points to implements log-Cholesky in my situation, but I am simply not able to see where is the Choleksy decomposition to be implemented in this my workflow. > > Regards, > B.Nataraj > > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Sunday, June 17, 2012 4:10 AM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see below. > > On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: >> Dear Mr.Kjetil, >> >> Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. >> >> I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. >> >> Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. > > I dont think there is a prewritten function for log cholesy, but it is > very easy to write! > >> makelogchol > function(A) {#A should be a square posdef matrix > m <- dim(A)[1] > uind <- upper.tri(A) > R <- chol(A) > diag <- log(diag(R)) > upper <- R[uind] > return(list(diag=diag, upper=upper)) > } > > and its "inverse" --- to get the cholesky factor back: > >> makemat > function(diag, upper) { > m <- length(diag) > mu <- m*(m-1)/2 > if (length(upper)!=mu) stop("incompatible lengths") > A <- matrix(0.0, m, m) > ind <- upper.tri(A) > A[ind] <- upper > diag(A) <- exp(diag) > A > } > > and an example of its use: > > Lets say A is the posdef matrix where you will start your > optimization, so you need > the log cholesy parameters to feed to optim: > >> A >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> start <- makelogchol(A) >> start > \$diag > [1] 0.5493061 0.4904146 0.4581454 > > \$upper > [1] 0.5773503 0.5773503 0.4082483 > > Then to get the cholesky factors back you call makemat: > >> R <- makemat(start\$diag,  start\$upper) >> R >         [,1]      [,2]      [,3] > [1,] 1.732051 0.5773503 0.5773503 > [2,] 0.000000 1.6329932 0.4082483 > [3,] 0.000000 0.0000000 1.5811388 >> t(R) %*% R >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> > > Kjetil > > > > > > > >> >> Regards, >> B.Nataraj >> >> >> >> -----Original Message----- >> From: Kjetil Halvorsen [mailto:[hidden email]] >> Sent: Friday, June 15, 2012 11:09 PM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email]; [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> see inline. >> >> On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >>> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. >> >> ¿What are you trying to optimize? If the objective function depends on >> a matrix argument >> which has to be a positive definite function, you must parametrize the >> matrix such that the matrix >> inside the optimizer always is positive definite. So if your positive >> definite matrix is A, then, for example, represent it as >> its cholesky decomposition A= L L^T where L is lower triangular with >> positive diagonal. Here the stricly >> upper diagonal part varies freely, but the diagonal not, so represent >> the diagonal as exp( l_i) >> where now the l_i varies freely. This is called the log-Cholesky >> parametrization. For other ideas along this lines, see the paper by >> Douglas Bates:  "Unconstrained Parameterizations for >> Variance-Covariance Matrices  " >> which you can find by googling. >> >> Kjetil >> >> >>  This I observed when I change the maximum iteration of the optim >> function set to 1 and upto iteration no. 3 it runs , it stuck at >> iteration 4 and above. >>> >>> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >>> >>> Regards, >>> B.Nataraj >>> >>> >>> -----Original Message----- >>> From: Bert Gunter [mailto:[hidden email]] >>> Sent: Friday, June 15, 2012 1:51 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> Follow the posting guide,please: I believe at this point we need >>> reproducible code and your data to provide you help. See ?dput to post >>> your matrix. >>> >>> -- Bert >>> >>> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>>> >>>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>>> >>>>> Error in chol.default(M_cov) : >>>>>  the leading minor of order 10 is not positive definite >>>> >>>> is surfaced during the function call by optim. >>>> >>>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>>> >>>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>>> >>>> Regards, >>>> B.Nataraj >>>> >>>> >>>> >>>> >>>> >>>> -----Original Message----- >>>> From: Bert Gunter [mailto:[hidden email]] >>>> Sent: Thursday, June 14, 2012 6:18 PM >>>> To: Nataraj B (ORLL-Biotech) >>>> Cc: [hidden email] >>>> Subject: Re: [R] Cholesky decomposition error >>>> >>>> Your matrix is not symmetric, positive definite. If you don't know >>>> what this means, you shouldn't be using chol() >>>> >>>> This may be because it isn't to begin with, or due to numerical error, >>>> it doesn't behave as one in the decomposition. My relative ignorance >>>> of numeric methods for linear algebra prevents me from saying more >>>> than that. >>>> >>>> -- Bert >>>> >>>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>>> Dear friends, >>>>> >>>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>>> >>>>> Error in chol.default(M_cov) : >>>>>  the leading minor of order 10 is not positive definite >>>>> >>>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>>> >>>>> Could you please help me to find where things are going wrong in my matrix? >>>>> >>>>> >>>>> Thanks and regards, >>>>> B.Natarj >>>>> >>>>> ______________________________________________ >>>>> [hidden email] mailing list >>>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>>> and provide commented, minimal, self-contained, reproducible code. >>>> >>>> >>>> >>>> -- >>>> >>>> Bert Gunter >>>> Genentech Nonclinical Biostatistics >>>> >>>> Internal Contact Info: >>>> Phone: 467-7374 >>>> Website: >>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>> >>>> >>> >>> >>> >>> -- >>> >>> Bert Gunter >>> Genentech Nonclinical Biostatistics >>> >>> Internal Contact Info: >>> Phone: 467-7374 >>> Website: >>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>> >>> ______________________________________________ >>> [hidden email] mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>> and provide commented, minimal, self-contained, reproducible code. >> >> > >         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 On Jun 19, 2012, at 10:31 , <[hidden email]> <[hidden email]> wrote: > Dear Kjetil, thanks for your time to detail the code. Sorry, I am still not able to grasp why you are creating random normal distribution for the covariance matrix and also why I should optimize the whole matrix value while I only need to optimize the two variables ?. > > This is my 15 x 15 covariance matrix >> M_cov > [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] > [1,]    4    3    3    2    2    1    2    2    3     2     2     2     3     3     1 > [2,]    3    5    3    2    2    1    2    2    2     2     2     2     2     2     1 > [3,]    3    3    3    2    2    1    2    2    2     2     2     2     2     2     1 > [4,]    2    2    2    3    2    1    2    3    2     1     1     1     1     1     1 > [5,]    2    2    2    2    2    1    2    2    2     1     1     1     1     1     1 > [6,]    1    1    1    1    1    1    1    1    1     1     0     1     1     0     0 > [7,]    2    2    2    2    2    1    3    2    2     1     1     1     1     1     1 > [8,]    2    2    2    3    2    1    2    4    2     1     1     1     1     1     1 > [9,]    3    2    2    2    2    1    2    2    3     1     1     1     2     2     1 > [10,]    2    2    2    1    1    1    1    1    1     2     1     2     2     1     0 > [11,]    2    2    2    1    1    0    1    1    1     1     2     1     1     2     1 > [12,]    2    2    2    1    1    1    1    1    1     2     1     3     2     1     0 > [13,]    3    2    2    1    1    1    1    1    2     2     1     2     3     2     0 > [14,]    3    2    2    1    1    0    1    1    2     1     2     1     2     3     1 > [15,]    1    1    1    1    1    0    1    1    1     0     1     0     0     1     1 > Notice that that matrix has rank 9 since eigenvalues no. 10 and up are effectively zero: > eigen(M_cov) \$values  [1]  2.539980e+01  5.456591e+00  3.694394e+00  2.757612e+00  1.746257e+00  [6]  1.385264e+00  6.232218e-01  5.853408e-01  3.515233e-01  1.480189e-16 [11]  9.492860e-17 -5.918493e-17 -1.142008e-16 -2.348974e-16 -7.429798e-16 \$vectors ... > > And the function to check the posdef condition is > >> isPosDef <- function(M) > { if ( all(M == t(M) ) ) {  # first test symmetric-ity >  if (  all(eigen(M)\$values > 0) ) {TRUE} >    else {FALSE} } #  else {FALSE}  # not symmetric > } > > And my call for optimization is > >> optim(c(rep(1,2)),obj,method="BFGS",hessian=TRUE,F1=F1,mu1=mu1,M_cov=M_cov) > > And I pass the vector F1 and mu1 used in my objective function along with M_cov (the covariance matrix) > > My objective function is > > >> obj<-function(theta,F1,mu1,M_cov) { > Vs<-theta[1] > Vn<-theta[2] > k<-length(F1) > for(i in 1:15) { > for(j in 1:15) { > > # to make the covariance matrix with the variables to optimize > if(i != j) > { > M_cov[i,j]<-M_cov[i,j]*Vs > } > else > M_cov[i,j]<-M_cov[i,j]*Vs+Vn > } > } > ...alias M_cov <- Vs * M_cov diag(M_cov) <- diag(M_cov) + Vn Right? Anyways, immediate observations are that (1) Vn had better be positive (2) Vs could be negative but 25.993*Vs + Vn must be positive. A little thought indicates that these are also sufficient conditions for M_cov to be pos.def. So you can replace the test for positive definiteness with bounds on your parameters. My instinct tells me that you might not be able to interpret the Vs < 0 case, in which case, the easy way is to just check that Vs and Vn are positive. People often ensure this by using a log-parametrization, Vs = exp(lVs), etc. -pd > # here I check the posdef condition of the modified covariance matrix in #earlier step > checkPD<-isPosDef(M_cov) > if(checkPD) > { > e<-t(F1-mu1)%*%solve(M_cov)%*%(F1-mu1) > logl<- -k/2*log(2*pi)-0.5*log(det(M_cov))-0.5*e > } > > # negative value since the function is for maximization > return(-logl) > } > > Is there anything missing in my approach of function optimization for variable estimation of a multivariate normal distribution? I looking forward your valuable comments. > > Regards, > B.Nataraj > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Monday, June 18, 2012 10:07 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see inline! > > On Mon, Jun 18, 2012 at 12:38 AM,  <[hidden email]> wrote: >> Thanks Kjetil for your detailed code for log-Cholesky but my situation is like optimization of the variables inside the matrix. > > But you must of course call my function makemat inside the call > to optim() or whatever function you use for optimization! > > > Here is an extended example, assuming the functions I have defined earlier: > > library(MASS) # for mvrnorm > library(mvtnorm) # for dmvnorm > > Let us simulate a multinorma sample with expectation zero (to simplify > the example) > and covariance matrix Sigma: > >> Sigma >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> > >> testdat <- mvrnorm(n=100, rep(0, 3), Sigma) > > We can calculate the usual maxlik estimator of Sigma: > >> (t(testdat) %*% testdat) /100 >          [,1]      [,2]     [,3] > [1,] 3.2424618 0.9670133 1.673897 > [2,] 0.9670133 2.9447654 1.289127 > [3,] 1.6738970 1.2891268 3.084106 > >> objective <- function(Sigma) { >        sum(apply(testdat, 1, FUN=function(x) dmvnorm(x, rep(0, 3), > Sigma, log=TRUE))) >    } >> objective(Sigma) > >> start > \$diag > [1] 0.5493061 0.4904146 0.4581454 > > \$upper > [1] 0.5773503 0.5773503 0.4082483 > >> startSigma <- c(start\$diag, start\$upper) >> startSigma > [1] 0.5493061 0.4904146 0.4581454 0.5773503 0.5773503 0.4082483 > >> objective(Sigma) > [1] -571.5945 >> objective(Sigma/3) > [1] -699.0552 >> objective(3*Sigma) > [1] -638.9688 >> > > > Finally: > > optim(startSigma, fn=function(parvec) { >        R <- makemat(parvec[1:3], parvec[4:6]) >        S <- t(R) %*% R >        print(S)   # remove this when works! >       obj <- -objective(S) >       obj >    }, control=list(REPORT=1)) > > # with method="BFGS" dopes not work well,  because takes too long > # steps! default Nelder-Mead is slow,  but works. > > > Giving the following output: > <  much removed ...> > \$par > [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 > > \$value > [1] 567.648 > > \$counts > function gradient >     501       NA > > \$convergence > [1] 1 > > \$message > NULL > >>> res <- .Last.value >> res\$par > [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 >> makemat(res\$par[1:3], res\$par[4:6]) >         [,1]     [,2]      [,3] > [1,] 1.788146 0.546956 0.9248746 > [2,] 0.000000 1.643628 0.4896758 > [3,] 0.000000 0.000000 1.4170710 >> R <-makemat(res\$par[1:3], res\$par[4:6]) >> t(R) %*% R >          [,1]      [,2]     [,3] > [1,] 3.1974676 0.9780374 1.653811 > [2,] 0.9780374 3.0006752 1.310711 > [3,] 1.6538112 1.3107108 3.103266 >> Sigma >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 > > > > ¿Easy? > > >> >> Let me change your own example matrix to explain the problem. >> >> A martrix is      [,1]         [,2]            [,3] >>               [1,]    3*Vs+Vn    1*Vs         1*Vs >>               [2,]    1*Vs     3*Vs+Vn    1*Vs >>               [3,]    1*Vs     1*Vs       3*Vs+Vn > > So make your own makemat function specific to this problem! > which you then use inside the call to optim. > > kjetil >> >> So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood function, which need the determinant and inverse of the matrix A to compute the value. >> >> Isn't it logically correct that I have to seed some initial values for the Vs and Vn and use it for predicting the value of the matrix and check for its positive definite condition, > > If you construct the matrix by above method, it will automatically be > posdef, so no need to check! > (except as check on programming ... ) > > if the condition fulfill, then I have to compute the determinant value > and inverse of the matrix in order to use them in the maximum > likelihood function and the optimization iteration to be carried out > further until parameters converges. >> >> I sincerely believe that you have points to implements log-Cholesky in my situation, but I am simply not able to see where is the Choleksy decomposition to be implemented in this my workflow. >> >> Regards, >> B.Nataraj >> >> >> >> >> -----Original Message----- >> From: Kjetil Halvorsen [mailto:[hidden email]] >> Sent: Sunday, June 17, 2012 4:10 AM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> see below. >> >> On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: >>> Dear Mr.Kjetil, >>> >>> Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. >>> >>> I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. >>> >>> Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. >> >> I dont think there is a prewritten function for log cholesy, but it is >> very easy to write! >> >>> makelogchol >> function(A) {#A should be a square posdef matrix >> m <- dim(A)[1] >> uind <- upper.tri(A) >> R <- chol(A) >> diag <- log(diag(R)) >> upper <- R[uind] >> return(list(diag=diag, upper=upper)) >> } >> >> and its "inverse" --- to get the cholesky factor back: >> >>> makemat >> function(diag, upper) { >> m <- length(diag) >> mu <- m*(m-1)/2 >> if (length(upper)!=mu) stop("incompatible lengths") >> A <- matrix(0.0, m, m) >> ind <- upper.tri(A) >> A[ind] <- upper >> diag(A) <- exp(diag) >> A >> } >> >> and an example of its use: >> >> Lets say A is the posdef matrix where you will start your >> optimization, so you need >> the log cholesy parameters to feed to optim: >> >>> A >>    [,1] [,2] [,3] >> [1,]    3    1    1 >> [2,]    1    3    1 >> [3,]    1    1    3 >>> start <- makelogchol(A) >>> start >> \$diag >> [1] 0.5493061 0.4904146 0.4581454 >> >> \$upper >> [1] 0.5773503 0.5773503 0.4082483 >> >> Then to get the cholesky factors back you call makemat: >> >>> R <- makemat(start\$diag,  start\$upper) >>> R >>        [,1]      [,2]      [,3] >> [1,] 1.732051 0.5773503 0.5773503 >> [2,] 0.000000 1.6329932 0.4082483 >> [3,] 0.000000 0.0000000 1.5811388 >>> t(R) %*% R >>    [,1] [,2] [,3] >> [1,]    3    1    1 >> [2,]    1    3    1 >> [3,]    1    1    3 >>> >> >> Kjetil >> >> >> >> >> >> >> >>> >>> Regards, >>> B.Nataraj >>> >>> >>> >>> -----Original Message----- >>> From: Kjetil Halvorsen [mailto:[hidden email]] >>> Sent: Friday, June 15, 2012 11:09 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email]; [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> see inline. >>> >>> On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >>>> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. >>> >>> ¿What are you trying to optimize? If the objective function depends on >>> a matrix argument >>> which has to be a positive definite function, you must parametrize the >>> matrix such that the matrix >>> inside the optimizer always is positive definite. So if your positive >>> definite matrix is A, then, for example, represent it as >>> its cholesky decomposition A= L L^T where L is lower triangular with >>> positive diagonal. Here the stricly >>> upper diagonal part varies freely, but the diagonal not, so represent >>> the diagonal as exp( l_i) >>> where now the l_i varies freely. This is called the log-Cholesky >>> parametrization. For other ideas along this lines, see the paper by >>> Douglas Bates:  "Unconstrained Parameterizations for >>> Variance-Covariance Matrices  " >>> which you can find by googling. >>> >>> Kjetil >>> >>> >>> This I observed when I change the maximum iteration of the optim >>> function set to 1 and upto iteration no. 3 it runs , it stuck at >>> iteration 4 and above. >>>> >>>> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >>>> >>>> Regards, >>>> B.Nataraj >>>> >>>> >>>> -----Original Message----- >>>> From: Bert Gunter [mailto:[hidden email]] >>>> Sent: Friday, June 15, 2012 1:51 PM >>>> To: Nataraj B (ORLL-Biotech) >>>> Cc: [hidden email] >>>> Subject: Re: [R] Cholesky decomposition error >>>> >>>> Follow the posting guide,please: I believe at this point we need >>>> reproducible code and your data to provide you help. See ?dput to post >>>> your matrix. >>>> >>>> -- Bert >>>> >>>> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>>>> >>>>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>>>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>>>> >>>>>> Error in chol.default(M_cov) : >>>>>> the leading minor of order 10 is not positive definite >>>>> >>>>> is surfaced during the function call by optim. >>>>> >>>>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>>>> >>>>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>>>> >>>>> Regards, >>>>> B.Nataraj >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> -----Original Message----- >>>>> From: Bert Gunter [mailto:[hidden email]] >>>>> Sent: Thursday, June 14, 2012 6:18 PM >>>>> To: Nataraj B (ORLL-Biotech) >>>>> Cc: [hidden email] >>>>> Subject: Re: [R] Cholesky decomposition error >>>>> >>>>> Your matrix is not symmetric, positive definite. If you don't know >>>>> what this means, you shouldn't be using chol() >>>>> >>>>> This may be because it isn't to begin with, or due to numerical error, >>>>> it doesn't behave as one in the decomposition. My relative ignorance >>>>> of numeric methods for linear algebra prevents me from saying more >>>>> than that. >>>>> >>>>> -- Bert >>>>> >>>>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>>>> Dear friends, >>>>>> >>>>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>>>> >>>>>> Error in chol.default(M_cov) : >>>>>> the leading minor of order 10 is not positive definite >>>>>> >>>>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>>>> >>>>>> Could you please help me to find where things are going wrong in my matrix? >>>>>> >>>>>> >>>>>> Thanks and regards, >>>>>> B.Natarj >>>>>> >>>>>> ______________________________________________ >>>>>> [hidden email] mailing list >>>>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>>>> and provide commented, minimal, self-contained, reproducible code. >>>>> >>>>> >>>>> >>>>> -- >>>>> >>>>> Bert Gunter >>>>> Genentech Nonclinical Biostatistics >>>>> >>>>> Internal Contact Info: >>>>> Phone: 467-7374 >>>>> Website: >>>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>>> >>>>> >>>> >>>> >>>> >>>> -- >>>> >>>> Bert Gunter >>>> Genentech Nonclinical Biostatistics >>>> >>>> Internal Contact Info: >>>> Phone: 467-7374 >>>> Website: >>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>> >>>> ______________________________________________ >>>> [hidden email] mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> >> >> > > > [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Reply | Threaded
Open this post in threaded view
|
Report Content as Inappropriate

Re: Cholesky decomposition error

 Dear Peter, Thanks for your reply, as per my last post; my current problem is not on positive definite related but to compare the objective function narrated by Kjetil for optimization of parameters Vs Objective function that I have posted in my last post. Though I appreciate the Kjetil approach that he uses only the covariance matrix information to optimize the parameters of the matrix but in my case I have other values to be used in the objective function to optimize my matrix parameters and I can not figure out how the Kjetil approach extendable for my specific problem and in this regard my inexperience in R as well as lack of knowledge of linear algebra is to be blamed but I have to get out of the problem earlier, is my current situation :-). Regards, B.Nataraj -----Original Message----- From: peter dalgaard [mailto:[hidden email]] Sent: Tuesday, June 19, 2012 8:07 PM To: Nataraj B (ORLL-Biotech) Cc: [hidden email]; [hidden email] Subject: Re: [R] Cholesky decomposition error On Jun 19, 2012, at 10:31 , <[hidden email]> <[hidden email]> wrote: > Dear Kjetil, thanks for your time to detail the code. Sorry, I am still not able to grasp why you are creating random normal distribution for the covariance matrix and also why I should optimize the whole matrix value while I only need to optimize the two variables ?. > > This is my 15 x 15 covariance matrix >> M_cov > [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] > [1,]    4    3    3    2    2    1    2    2    3     2     2     2     3     3     1 > [2,]    3    5    3    2    2    1    2    2    2     2     2     2     2     2     1 > [3,]    3    3    3    2    2    1    2    2    2     2     2     2     2     2     1 > [4,]    2    2    2    3    2    1    2    3    2     1     1     1     1     1     1 > [5,]    2    2    2    2    2    1    2    2    2     1     1     1     1     1     1 > [6,]    1    1    1    1    1    1    1    1    1     1     0     1     1     0     0 > [7,]    2    2    2    2    2    1    3    2    2     1     1     1     1     1     1 > [8,]    2    2    2    3    2    1    2    4    2     1     1     1     1     1     1 > [9,]    3    2    2    2    2    1    2    2    3     1     1     1     2     2     1 > [10,]    2    2    2    1    1    1    1    1    1     2     1     2     2     1     0 > [11,]    2    2    2    1    1    0    1    1    1     1     2     1     1     2     1 > [12,]    2    2    2    1    1    1    1    1    1     2     1     3     2     1     0 > [13,]    3    2    2    1    1    1    1    1    2     2     1     2     3     2     0 > [14,]    3    2    2    1    1    0    1    1    2     1     2     1     2     3     1 > [15,]    1    1    1    1    1    0    1    1    1     0     1     0     0     1     1 > Notice that that matrix has rank 9 since eigenvalues no. 10 and up are effectively zero: > eigen(M_cov) \$values  [1]  2.539980e+01  5.456591e+00  3.694394e+00  2.757612e+00  1.746257e+00  [6]  1.385264e+00  6.232218e-01  5.853408e-01  3.515233e-01  1.480189e-16 [11]  9.492860e-17 -5.918493e-17 -1.142008e-16 -2.348974e-16 -7.429798e-16 \$vectors ... > > And the function to check the posdef condition is > >> isPosDef <- function(M) > { if ( all(M == t(M) ) ) {  # first test symmetric-ity >  if (  all(eigen(M)\$values > 0) ) {TRUE} >    else {FALSE} } #  else {FALSE}  # not symmetric > } > > And my call for optimization is > >> optim(c(rep(1,2)),obj,method="BFGS",hessian=TRUE,F1=F1,mu1=mu1,M_cov=M_cov) > > And I pass the vector F1 and mu1 used in my objective function along with M_cov (the covariance matrix) > > My objective function is > > >> obj<-function(theta,F1,mu1,M_cov) { > Vs<-theta[1] > Vn<-theta[2] > k<-length(F1) > for(i in 1:15) { > for(j in 1:15) { > > # to make the covariance matrix with the variables to optimize > if(i != j) > { > M_cov[i,j]<-M_cov[i,j]*Vs > } > else > M_cov[i,j]<-M_cov[i,j]*Vs+Vn > } > } > ...alias M_cov <- Vs * M_cov diag(M_cov) <- diag(M_cov) + Vn Right? Anyways, immediate observations are that (1) Vn had better be positive (2) Vs could be negative but 25.993*Vs + Vn must be positive. A little thought indicates that these are also sufficient conditions for M_cov to be pos.def. So you can replace the test for positive definiteness with bounds on your parameters. My instinct tells me that you might not be able to interpret the Vs < 0 case, in which case, the easy way is to just check that Vs and Vn are positive. People often ensure this by using a log-parametrization, Vs = exp(lVs), etc. -pd > # here I check the posdef condition of the modified covariance matrix in #earlier step > checkPD<-isPosDef(M_cov) > if(checkPD) > { > e<-t(F1-mu1)%*%solve(M_cov)%*%(F1-mu1) > logl<- -k/2*log(2*pi)-0.5*log(det(M_cov))-0.5*e > } > > # negative value since the function is for maximization > return(-logl) > } > > Is there anything missing in my approach of function optimization for variable estimation of a multivariate normal distribution? I looking forward your valuable comments. > > Regards, > B.Nataraj > > > > -----Original Message----- > From: Kjetil Halvorsen [mailto:[hidden email]] > Sent: Monday, June 18, 2012 10:07 PM > To: Nataraj B (ORLL-Biotech) > Cc: [hidden email] > Subject: Re: [R] Cholesky decomposition error > > see inline! > > On Mon, Jun 18, 2012 at 12:38 AM,  <[hidden email]> wrote: >> Thanks Kjetil for your detailed code for log-Cholesky but my situation is like optimization of the variables inside the matrix. > > But you must of course call my function makemat inside the call > to optim() or whatever function you use for optimization! > > > Here is an extended example, assuming the functions I have defined earlier: > > library(MASS) # for mvrnorm > library(mvtnorm) # for dmvnorm > > Let us simulate a multinorma sample with expectation zero (to simplify > the example) > and covariance matrix Sigma: > >> Sigma >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 >> > >> testdat <- mvrnorm(n=100, rep(0, 3), Sigma) > > We can calculate the usual maxlik estimator of Sigma: > >> (t(testdat) %*% testdat) /100 >          [,1]      [,2]     [,3] > [1,] 3.2424618 0.9670133 1.673897 > [2,] 0.9670133 2.9447654 1.289127 > [3,] 1.6738970 1.2891268 3.084106 > >> objective <- function(Sigma) { >        sum(apply(testdat, 1, FUN=function(x) dmvnorm(x, rep(0, 3), > Sigma, log=TRUE))) >    } >> objective(Sigma) > >> start > \$diag > [1] 0.5493061 0.4904146 0.4581454 > > \$upper > [1] 0.5773503 0.5773503 0.4082483 > >> startSigma <- c(start\$diag, start\$upper) >> startSigma > [1] 0.5493061 0.4904146 0.4581454 0.5773503 0.5773503 0.4082483 > >> objective(Sigma) > [1] -571.5945 >> objective(Sigma/3) > [1] -699.0552 >> objective(3*Sigma) > [1] -638.9688 >> > > > Finally: > > optim(startSigma, fn=function(parvec) { >        R <- makemat(parvec[1:3], parvec[4:6]) >        S <- t(R) %*% R >        print(S)   # remove this when works! >       obj <- -objective(S) >       obj >    }, control=list(REPORT=1)) > > # with method="BFGS" dopes not work well,  because takes too long > # steps! default Nelder-Mead is slow,  but works. > > > Giving the following output: > <  much removed ...> > \$par > [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 > > \$value > [1] 567.648 > > \$counts > function gradient >     501       NA > > \$convergence > [1] 1 > > \$message > NULL > >>> res <- .Last.value >> res\$par > [1] 0.5811796 0.4969062 0.3485920 0.5469560 0.9248746 0.4896758 >> makemat(res\$par[1:3], res\$par[4:6]) >         [,1]     [,2]      [,3] > [1,] 1.788146 0.546956 0.9248746 > [2,] 0.000000 1.643628 0.4896758 > [3,] 0.000000 0.000000 1.4170710 >> R <-makemat(res\$par[1:3], res\$par[4:6]) >> t(R) %*% R >          [,1]      [,2]     [,3] > [1,] 3.1974676 0.9780374 1.653811 > [2,] 0.9780374 3.0006752 1.310711 > [3,] 1.6538112 1.3107108 3.103266 >> Sigma >     [,1] [,2] [,3] > [1,]    3    1    1 > [2,]    1    3    1 > [3,]    1    1    3 > > > > ¿Easy? > > >> >> Let me change your own example matrix to explain the problem. >> >> A martrix is      [,1]         [,2]            [,3] >>               [1,]    3*Vs+Vn    1*Vs         1*Vs >>               [2,]    1*Vs     3*Vs+Vn    1*Vs >>               [3,]    1*Vs     1*Vs       3*Vs+Vn > > So make your own makemat function specific to this problem! > which you then use inside the call to optim. > > kjetil >> >> So I have to optimize two variables "Vs" and "Vn" using the maximum likelihood function, which need the determinant and inverse of the matrix A to compute the value. >> >> Isn't it logically correct that I have to seed some initial values for the Vs and Vn and use it for predicting the value of the matrix and check for its positive definite condition, > > If you construct the matrix by above method, it will automatically be > posdef, so no need to check! > (except as check on programming ... ) > > if the condition fulfill, then I have to compute the determinant value > and inverse of the matrix in order to use them in the maximum > likelihood function and the optimization iteration to be carried out > further until parameters converges. >> >> I sincerely believe that you have points to implements log-Cholesky in my situation, but I am simply not able to see where is the Choleksy decomposition to be implemented in this my workflow. >> >> Regards, >> B.Nataraj >> >> >> >> >> -----Original Message----- >> From: Kjetil Halvorsen [mailto:[hidden email]] >> Sent: Sunday, June 17, 2012 4:10 AM >> To: Nataraj B (ORLL-Biotech) >> Cc: [hidden email] >> Subject: Re: [R] Cholesky decomposition error >> >> see below. >> >> On Fri, Jun 15, 2012 at 11:53 PM,  <[hidden email]> wrote: >>> Dear Mr.Kjetil, >>> >>> Thanks for your comment. You have already pointed me the article in reply to one of my earlier post to this list and I am following the paper. Now I am checking for condition for positive definiteness for original matrix using a simple script (got from earlier posting in the list) and if the condition passed then the optimization function perform decomposition of the matrix using chol() function. >>> >>> I believe all 5 parameterization techniques listed in the paper weigh each other based on its performance and not on accuracy, please correct me if I am wrong. Since I am novice, I just want to start with the easiest method "Cholesky" for my purpose. >>> >>> Are you recommending log-Cholesky function , if so could you please tell me the name of the function which does log-Cholesky decomposition of a matrix. >> >> I dont think there is a prewritten function for log cholesy, but it is >> very easy to write! >> >>> makelogchol >> function(A) {#A should be a square posdef matrix >> m <- dim(A)[1] >> uind <- upper.tri(A) >> R <- chol(A) >> diag <- log(diag(R)) >> upper <- R[uind] >> return(list(diag=diag, upper=upper)) >> } >> >> and its "inverse" --- to get the cholesky factor back: >> >>> makemat >> function(diag, upper) { >> m <- length(diag) >> mu <- m*(m-1)/2 >> if (length(upper)!=mu) stop("incompatible lengths") >> A <- matrix(0.0, m, m) >> ind <- upper.tri(A) >> A[ind] <- upper >> diag(A) <- exp(diag) >> A >> } >> >> and an example of its use: >> >> Lets say A is the posdef matrix where you will start your >> optimization, so you need >> the log cholesy parameters to feed to optim: >> >>> A >>    [,1] [,2] [,3] >> [1,]    3    1    1 >> [2,]    1    3    1 >> [3,]    1    1    3 >>> start <- makelogchol(A) >>> start >> \$diag >> [1] 0.5493061 0.4904146 0.4581454 >> >> \$upper >> [1] 0.5773503 0.5773503 0.4082483 >> >> Then to get the cholesky factors back you call makemat: >> >>> R <- makemat(start\$diag,  start\$upper) >>> R >>        [,1]      [,2]      [,3] >> [1,] 1.732051 0.5773503 0.5773503 >> [2,] 0.000000 1.6329932 0.4082483 >> [3,] 0.000000 0.0000000 1.5811388 >>> t(R) %*% R >>    [,1] [,2] [,3] >> [1,]    3    1    1 >> [2,]    1    3    1 >> [3,]    1    1    3 >>> >> >> Kjetil >> >> >> >> >> >> >> >>> >>> Regards, >>> B.Nataraj >>> >>> >>> >>> -----Original Message----- >>> From: Kjetil Halvorsen [mailto:[hidden email]] >>> Sent: Friday, June 15, 2012 11:09 PM >>> To: Nataraj B (ORLL-Biotech) >>> Cc: [hidden email]; [hidden email] >>> Subject: Re: [R] Cholesky decomposition error >>> >>> see inline. >>> >>> On Fri, Jun 15, 2012 at 4:33 AM,  <[hidden email]> wrote: >>>> Thanks for your reply. I am sorry and I am bit hurried up to say before doing a proper due diligence, I have found out that during the optimization the variables tend to vary the values of the matrix , the function report error at some point (in particular iteration step) when the matrix become non-decomposable due to not a positive definiteness. >>> >>> ¿What are you trying to optimize? If the objective function depends on >>> a matrix argument >>> which has to be a positive definite function, you must parametrize the >>> matrix such that the matrix >>> inside the optimizer always is positive definite. So if your positive >>> definite matrix is A, then, for example, represent it as >>> its cholesky decomposition A= L L^T where L is lower triangular with >>> positive diagonal. Here the stricly >>> upper diagonal part varies freely, but the diagonal not, so represent >>> the diagonal as exp( l_i) >>> where now the l_i varies freely. This is called the log-Cholesky >>> parametrization. For other ideas along this lines, see the paper by >>> Douglas Bates:  "Unconstrained Parameterizations for >>> Variance-Covariance Matrices  " >>> which you can find by googling. >>> >>> Kjetil >>> >>> >>> This I observed when I change the maximum iteration of the optim >>> function set to 1 and upto iteration no. 3 it runs , it stuck at >>> iteration 4 and above. >>>> >>>> Now, I am trying to find ways to escalate such a condition inside the function during the iteration process and if possible please help me to do that. >>>> >>>> Regards, >>>> B.Nataraj >>>> >>>> >>>> -----Original Message----- >>>> From: Bert Gunter [mailto:[hidden email]] >>>> Sent: Friday, June 15, 2012 1:51 PM >>>> To: Nataraj B (ORLL-Biotech) >>>> Cc: [hidden email] >>>> Subject: Re: [R] Cholesky decomposition error >>>> >>>> Follow the posting guide,please: I believe at this point we need >>>> reproducible code and your data to provide you help. See ?dput to post >>>> your matrix. >>>> >>>> -- Bert >>>> >>>> On Thu, Jun 14, 2012 at 11:30 PM,  <[hidden email]> wrote: >>>>> >>>>> Thanks for your reply. To my surprise I can find one more strange behavior of  my 15X15 matrix "A", that is if I call the function  chol(A) in the terminal it decompose the matrix fine without any errors or warnings. >>>>> But if I call the function chol() within a function, which I have written in order to call the function (contains formula) for optimization routine "optim()" and also supplied with the same matrix "A" as argument, the error mentioned >>>>> >>>>>> Error in chol.default(M_cov) : >>>>>> the leading minor of order 10 is not positive definite >>>>> >>>>> is surfaced during the function call by optim. >>>>> >>>>> Why the matrix fulfill the symmetric and positive definite for chol() in one case but fails in other case when the function chol() is called in other function ? >>>>> >>>>> I played around parameters of "optim" function but nothing seems to be working and I am confused and I am looking for some hints to introspect the problem further. >>>>> >>>>> Regards, >>>>> B.Nataraj >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> -----Original Message----- >>>>> From: Bert Gunter [mailto:[hidden email]] >>>>> Sent: Thursday, June 14, 2012 6:18 PM >>>>> To: Nataraj B (ORLL-Biotech) >>>>> Cc: [hidden email] >>>>> Subject: Re: [R] Cholesky decomposition error >>>>> >>>>> Your matrix is not symmetric, positive definite. If you don't know >>>>> what this means, you shouldn't be using chol() >>>>> >>>>> This may be because it isn't to begin with, or due to numerical error, >>>>> it doesn't behave as one in the decomposition. My relative ignorance >>>>> of numeric methods for linear algebra prevents me from saying more >>>>> than that. >>>>> >>>>> -- Bert >>>>> >>>>> On Thu, Jun 14, 2012 at 4:23 AM,  <[hidden email]> wrote: >>>>>> Dear friends, >>>>>> >>>>>> When I do Cholesky decomposition for a 15x15 matrix using the function chol(), I get the following error for which I do not understand the meaning of the error >>>>>> >>>>>> Error in chol.default(M_cov) : >>>>>> the leading minor of order 10 is not positive definite >>>>>> >>>>>> When I searched online for similar error reported earlier I could get few hits but not of much help to resolve my error and one post suggested to use different function called sechol() from accuracy package but that did not work and it leads to different errors. So I want to stick to function chol() itself. >>>>>> >>>>>> Could you please help me to find where things are going wrong in my matrix? >>>>>> >>>>>> >>>>>> Thanks and regards, >>>>>> B.Natarj >>>>>> >>>>>> ______________________________________________ >>>>>> [hidden email] mailing list >>>>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>>>> and provide commented, minimal, self-contained, reproducible code. >>>>> >>>>> >>>>> >>>>> -- >>>>> >>>>> Bert Gunter >>>>> Genentech Nonclinical Biostatistics >>>>> >>>>> Internal Contact Info: >>>>> Phone: 467-7374 >>>>> Website: >>>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>>> >>>>> >>>> >>>> >>>> >>>> -- >>>> >>>> Bert Gunter >>>> Genentech Nonclinical Biostatistics >>>> >>>> Internal Contact Info: >>>> Phone: 467-7374 >>>> Website: >>>> http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm>>>> >>>> ______________________________________________ >>>> [hidden email] mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-help>>>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html>>>> and provide commented, minimal, self-contained, reproducible code. >>> >>> >> >> > > >       [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. -- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
12
Loading...