[R] Positive Definite Matrix

[David Winsemius]

On Jan 29, 2011, at 10:11 AM, David Winsemius wrote:

>
> 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.

In addition to the two "is." functions cited earlier there is also a  
"posdefify" function by Maechler in the sfsmisc package: "  
Description : From a matrix m, construct a "close" positive definite  
one."

-- 
David.
>>
>>>
>>>
>>> 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?
>>>>