[R] Positive Definite Matrix
Prof Brian Ripley
ripley at stats.ox.ac.uk
Sat Jan 29 22:00:19 CET 2011
On Sat, 29 Jan 2011, David Winsemius wrote:
>
> On Jan 29, 2011, at 12:17 PM, Prof Brian Ripley wrote:
>
>> On Sat, 29 Jan 2011, David Winsemius wrote:
>>
>>>
>>> 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."
>>
>> But again, that is not usually what you want. There is no guarantee that
>> the result is positive-definite enough that the Cholesky decomposition will
>> work.
>
> I don't see a Cholesky decomposition method being used in that function. It
> appears to my reading to be following what would be called an
> eigendecomposition.
Correct, but my point is that one does not usually want a
"close" positive definite one
but a 'square root'.
>
> --
> David.
>
>
>> Give up on Cholesky factors unless you have a matrix you know must be
>> symmetric and strictly positive definite, and use the eigendecomposition
>> instead (setting negative eigenvalues to zero). You can then work with the
>> factorization to ensure that (for example) variances are always
>> non-negative because they are always computed as sums of squares.
>>
>> This sort of thing is done in many of the multivariate analysis
>> calculations in R (e.g. cmdscale) and in well-designed packages.
>>
>>>
>>> --
>>> 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?
>>>
>>> ______________________________________________
>>> 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.
>>
>> --
>> Brian D. Ripley, ripley at stats.ox.ac.uk
>> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
>> University of Oxford, Tel: +44 1865 272861 (self)
>> 1 South Parks Road, +44 1865 272866 (PA)
>> Oxford OX1 3TG, UK Fax: +44 1865 272595
>
> David Winsemius, MD
> West Hartford, CT
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list