# [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.
>>>>>>> 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)
>>>>>>> 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
>> 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

```