[R] Positive Definite Matrix

[David Winsemius]

On Jan 29, 2011, at 9:59 AM, John Fox wrote:

> Dear David and Alex,
>
> I'd be a little careful about testing exact equality as in all(M ==  
> t(M) and
> careful as well about a test such as all(eigen(M)$values > 0) since  
> real
> arithmetic on a computer can't be counted on to be exact.

Which was why I pointed to that thread from 2005 and the existing work  
that had been put into packages. If you want to substitute all.equal  
for all, there might be fewer numerical false alarms, but I would  
think there could be other potential problems that might deserve  
warnings.


>
> Best,
> John
>
> --------------------------------
> John Fox
> Senator William McMaster
>  Professor of Social Statistics
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.mcmaster.ca/jfox
>
>
>
>
>> -----Original Message-----
>> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org 
>> ]
>> On Behalf Of David Winsemius
>> Sent: January-29-11 9:46 AM
>> To: Alex Smith
>> Cc: r-help at r-project.org Help
>> Subject: Re: [R] Positive Definite Matrix
>>
>>
>> On Jan 29, 2011, at 7:58 AM, David Winsemius wrote:
>>
>>>
>>> On Jan 29, 2011, at 7:22 AM, Alex Smith wrote:
>>>
>>>> Hello I am trying to determine wether a given matrix is symmetric  
>>>> and
>>>> positive matrix. The matrix has real valued elements.
>>>>
>>>> I have been reading about the cholesky method and another method is
>>>> to find the eigenvalues. I cant understand how to implement  
>>>> either of
>>>> the two. Can someone point me to the right direction. I have used
>>>> ?chol to see the help but if the matrix is not positive definite it
>>>> comes up as error. I know how to the get the eigenvalues but how  
>>>> can
>>>> I then put this into a program to check them as the just come up  
>>>> with
>>>> $values.
>>>>
>>>> Is checking that the eigenvalues are positive enough to determine
>>>> wether the matrix is positive definite?
>>>
>>> That is a fairly simple linear algebra fact that googling or pulling
>>> out a standard reference should have confirmed.
>>
>> Just to be clear (since on the basis of some off-line  
>> communications it
>> did not seem to be clear):  A real, symmetric matrix is Hermitian  
>> (and
>> therefore all of its eigenvalues are real). Further, it is positive-
>> definite if and only if its eigenvalues are all positive.
>>
>> qwe<-c(2,-1,0,-1,2,-1,0,1,2)
>> q<-matrix(qwe,nrow=3)
>>
>> 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
>>
>>                           }
>>
>>> isPosDef(q)
>> [1] FALSE
>>
>>>
>>>>
>>>> m
>>>>   [,1] [,2] [,3] [,4] [,5]
>>>> [1,]  1.0  0.0  0.5 -0.3  0.2
>>>> [2,]  0.0  1.0  0.1  0.0  0.0
>>>> [3,]  0.5  0.1  1.0  0.3  0.7
>>>> [4,] -0.3  0.0  0.3  1.0  0.4
>>>> [5,]  0.2  0.0  0.7  0.4  1.0
>>
>>> isPosDef(m)
>> [1] TRUE
>>
>> You might want to look at prior postings by people more  
>> knowledgeable than
>> me:
>>
>> http://finzi.psych.upenn.edu/R/Rhelp02/archive/57794.html
>>
>> Or look at what are probably better solutions in available packages:
>>
>> http://finzi.psych.upenn.edu/R/library/corpcor/html/rank.condition.html
>> http://finzi.psych.upenn.edu/R/library/matrixcalc/html/is.positive.definit
>> e.html
>>
>>
>> --
>> David.
>>
>>>>
>>>> this is the matrix that I know is positive definite.
>>>>
>>>> eigen(m)
>>>> $values
>>>> [1] 2.0654025 1.3391291 1.0027378 0.3956079 0.1971228
>>>>
>>>> $vectors
>>>>          [,1]        [,2]         [,3]        [,4]        [,5]
>>>> [1,] -0.32843233  0.69840166  0.080549876  0.44379474  0.44824689
>>>> [2,] -0.06080335  0.03564769 -0.993062427 -0.01474690  0.09296096
>>>> [3,] -0.64780034  0.12089168 -0.027187620  0.08912912 -0.74636235
>>>> [4,] -0.31765040 -0.68827876  0.007856812  0.60775962  0.23651023
>>>> [5,] -0.60653780 -0.15040584  0.080856897 -0.65231358  0.42123526
>>>>
>>>> and this are the eigenvalues and eigenvectors.
>>>> I thought of using
>>>> eigen(m,only.values=T)
>>>> $values
>>>> [1] 2.0654025 1.3391291 1.0027378 0.3956079 0.1971228
>>>>
>>>> $vectors
>>>> NULL
>>>>
>>>
>>>> m <- matrix(scan(textConnection("
>>>  1.0  0.0  0.5 -0.3  0.2
>>>  0.0  1.0  0.1  0.0  0.0
>>>  0.5  0.1  1.0  0.3  0.7
>>> -0.3  0.0  0.3  1.0  0.4
>>>  0.2  0.0  0.7  0.4  1.0
>>> ")), 5, byrow=TRUE)
>>> #Read 25 items
>>>> m
>>>    [,1] [,2] [,3] [,4] [,5]
>>> [1,]  1.0  0.0  0.5 -0.3  0.2
>>> [2,]  0.0  1.0  0.1  0.0  0.0
>>> [3,]  0.5  0.1  1.0  0.3  0.7
>>> [4,] -0.3  0.0  0.3  1.0  0.4
>>> [5,]  0.2  0.0  0.7  0.4  1.0
>>>
>>> all( eigen(m)$values >0 )
>>> #[1] TRUE
>>>
>>>> Then i thought of using logical expression to determine if there  
>>>> are
>>>> negative eigenvalues but couldnt work. I dont know what error  
>>>> this is
>>>>
>>>> b<-(a<0)
>>>> Error: (list) object cannot be coerced to type 'double'
>>>
>>> ??? where did "a" and "b" come from?
>>>
>>>>
>>>
>>>
>>> David Winsemius, MD
>>> West Hartford, CT
>>>
>>> ______________________________________________
>>> 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.
>>
>> David Winsemius, MD
>> West Hartford, CT
>>
>> ______________________________________________
>> 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.
>

David Winsemius, MD
West Hartford, CT