[R] Positive Definite Matrix

[John Fox]

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