[R] Cholesky decomposition error

Kjetil Halvorsen kjetilbrinchmannhalvorsen at gmail.com
Sun Jun 17 00:39:31 CEST 2012


see below.

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



More information about the R-help mailing list