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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] > On Behalf Of David Winsemius > Sent: January-29-11 9:46 AM > To: Alex Smith > Cc: r-help@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@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@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.
