see below.
On Fri, Jun 15, 2012 at 11:53 PM, <nata...@orchidpharma.com> 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:kjetilbrinchmannhalvor...@gmail.com]
> Sent: Friday, June 15, 2012 11:09 PM
> To: Nataraj B (ORLL-Biotech)
> Cc: gunter.ber...@gene.com; r-help@r-project.org
> Subject: Re: [R] Cholesky decomposition error
>
> see inline.
>
> On Fri, Jun 15, 2012 at 4:33 AM, <nata...@orchidpharma.com> 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:gunter.ber...@gene.com]
>> Sent: Friday, June 15, 2012 1:51 PM
>> To: Nataraj B (ORLL-Biotech)
>> Cc: r-help@r-project.org
>> 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, <nata...@orchidpharma.com> 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:gunter.ber...@gene.com]
>>> Sent: Thursday, June 14, 2012 6:18 PM
>>> To: Nataraj B (ORLL-Biotech)
>>> Cc: r-help@r-project.org
>>> 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, <nata...@orchidpharma.com> 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
>>>>
>>>> ______________________________________________
>>>> R-help@r-project.org 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
>>
>> ______________________________________________
>> R-help@r-project.org 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.
>
>
______________________________________________
R-help@r-project.org 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.
