Positive Definite Matrix

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.

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

David Winsemius, MD
West Hartford, CT