[R] Positive Definite Matrix
David Winsemius
dwinsemius at comcast.net
Sat Jan 29 20:02:27 CET 2011
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.
--
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
More information about the R-help
mailing list